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1:  (a)  Map  shows  the  epicenter  (black  star)  of  a  hypothetical  seismic  event,  whose 
synthetic  waveforms  at  Berkeley  Digital  Seismic  Network  stations  (black  triangles)  were 
used  for  evaluation  of  the  inversion  procedure  proposed  in  this  study,  assuming  different 
source-types,  (b)  Three-component  synthetic  displacement  wavefonns  (0.02-0.05  Hz)  for 
one  example  random  MT  solution  assuming  a  DC  source-type.  R  =  epicentral  distance, 

Az  =  azimuth,  and  Dmax  =  maximum  displacement  amplitude  at  a  station.  Beach-ball 
represents  lower  hemisphere  P-wave  radiation  pattern . 6 

2:  Best- fitting  MT  solutions  of  event  TE1  described  by  normalized  eigenvalues  (A), 
eigenvector  parameters  (a1,a2,a3,b1,b2)  and  seismic  moment  scale  factor  (ml)  for 
various  source-types:  (a)  tensile  crack  in  sediment,  (b)  tensile  crack  in  salt,  (c)  pure 
explosion  and  (d)  a  DC.  The  solid  and  dashed  lines  are  observed  and  synthetic 
displacement  waveforms,  respectively.  Meaning  of  other  symbols  is  same  as  in  Figure 
lb.  Station-specific  R,  Az  and  Dmax  are  same  in  all  subplots.  Waveforms  were  filtered  in 
the  pass  band  0. 1-0.3  Hz  for  stations  LA01,  LA02,  LA03  and  LA09,  and  in  0. 1-0.2  Hz 
for  station  LA08.  Final  values  of  (at,  a2,  a3,  blt  b2)  are  non-unique  (see  Figure  A2, 
Appendix) . 8 

3:  Grid  of  7457  unique  normalized  eigenvalues  (black  '+’  signs)  or  unique  source-types 
on  (a)  the  fundamental  Lune  (Tape  and  Tape,  20 12a, b),  and  (b)  on  the  Hudson  source- 
type  plot  (Hudson  et  al.,  1989).  Black  crosses  are  positions  of  major  theoretical  source- 
types  shown  with  their  un-nonnalized  eigenvalues . 10 

4:  NSS  of  the  three  events  using  low-frequency  displacement  wavefonns:  (a)  TE1,  (b) 

TE2  and  (c)  TE3.  Left  panels  (NSS  Inversion)  show  NSS  computed  using  the  inversion 
approach  in  this  study.  Right  panels  (NSS  Random)  show  NSS  computed  using  randomly 
generated  80  million  MT.  The  contours  and  shading  represent  absolute  values  of  VR  (%) 
whereas  NSS  plots  in  other  studies  usually  show  nonnalized  VR  (e.g.,  Guilhem  et  al, 

2014;  Nayak  and  Dreger,  2014;  Chiang  et  al.,  2014;  Boyd  et  al.,  2015).  For  each  event, 
the  VR  scale  is  same  for  both  NSS  plots  (left  and  right)  to  enable  better  comparison. 

Black  crosses  are  positions  of  major  theoretical  source-types.  For  each  event  the  white 
star  is  the  position  of  the  best-fitting  full  MT  solution  from  a  time-domain  full  MT 
inversion  of  wavefonns.  In  each  plot  the  white  circle  is  the  source-type  corresponding  to 
the  maximum  VR  recovered  by  each  NSS  (VRmax  in  the  lower  left  comer).  The  color 


version  of  this  figure  is  available  only  in  the  electronic  edition . 12 

5:  Figure  showing  the  true  sign  function  sign(x)  (equation  16),  its  approximation  used  in 
this  study,  uFM(x)  (equation  17),  and  its  normalized  derivative  dliF^(x\  where  x  =  itoz . 15 
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6:  NSS  of  event  TE3  using  P-wave  first  motion  polarities.  ‘+Dipole’  and  ‘+Crack’  are 
abbreviated  to  ‘+D’  and  ‘+C’,  respectively  in  the  left  panel.  Beach-balls  represent  P-wave 
lower  hemisphere  radiation  pattern  predicted  by  MT  solution  corresponding  to  the 
maximum  VR  recovered  by  each  NSS  (white  circle).  In  the  beach-balls,  black  crosses 
and  circles  represent  observed  positive  and  negative  P-wave  FM  polarities,  respectively; 
the  size  of  the  polarity  symbols  is  scaled  by  their  quality  weight  (1,  2  or  3).  ‘P’  and  ‘T’ 
indicate  pressure  and  Tension  axes,  respectively.  Explanation  of  other  symbols,  shading 
and  contours  in  NSS  plots  is  same  as  in  Figure  4.  The  color  version  of  this  figure  is 
available  only  in  the  electronic  edition .  17 

7:  NSS  of  the  three  events  using  both  low-frequency  displacement  waveforms  and  FM 
polarity  data:  (a)  TE1,  (b)  TE2  and  (c)  TE3.  Explanation  of  symbols,  shading  and 
contours  is  same  as  in  Figure  4.  The  color  version  of  this  figure  is  available  only  in  the 
electronic  edition .  19 

8:  a)  Green’s  functions  (GF)  computed  using  the  Song  et  ah,  (1996)  ID  western  U.S. 
velocity  model.  GF  are  in  displacement  and  filtered  between  10-50  second  period.  The 
traces  for  each  component,  from  top  to  bottom,  are  at  200,  400,  600,  800,  1000,  and  1200 
m  source  depths,  b)  Averaged  spectral  displacement  amplitudes  between  10  to  50 
seconds  for  the  ten  fundamental  Green’s  functions,  c)  Moment  tensor  elements  computed 


using  the  same  velocity  model  and  filtered  between  10-50  second  period . 25 

9:  Velocity  models  derived  from  the  Song  et  al.  (1996)  ID  model  by  keeping  the  top  2.5- 
km  vertical  travel  time  constant.  The  model  parameters  are  the  same  below  2.5-km  depth . 27 


10:  Isotropic  moment  and  total  seismic  moment  percent  change  for  a)  a  pure  explosion 
source  and  b)  a  composite  source,  plotted  as  a  function  of  source  depth.  The  average 
value  of  all  59  models  is  the  solid  red  line,  the  shaded  region  is  2-sigma  from  the  mean, 
and  the  dashed  line  represents  no  deviation  from  the  input  seismic  moments . 28 

11:  Fixed  Source  Depth  Sensitivity  Analysis.  Source  depth  is  fixed  while  we  tested  the 
full  suite  of  Green’s  function  depths  (x-axis),  and  compare  the  VR,  total  moment  percent 
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depth  for  (a)  a  pure  explosion  source  and  (b)  a  composite  source.  Similarly,  the  solid 
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13:  HUMMING  ALBATROSS  moment  tensor  solutions  from  time  domain  full  and 
deviatoric  inversion,  and  best  solution  from  combined  waveform  and  P-wave  first  motion 


(filled  circle)  Network  Sensitivity  Solutions.  All  polarities  are  up . 31 

14:  Production  shot  focal  mechanisms  and  waveform  fits  from  full  and  deviatoric 
moment  tensor  inversions . 32 


15:  Network  Sensitivity  Solutions  (NSS)  using  a)  waveforms,  b)  waveform  and  first 
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16:  Production  shot  moment  tensor  solutions  as  a  function  source  depth,  using  only 
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17:  Production  shot  full  moment  tensor  solutions  as  a  function  of  filter  passband,  using 
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18:  a)  Isotropic  moment  vs.  yield  from  nuclear  explosions  listed  in  Stevens  and  Murphy 
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inversion  vs.  the  true  yield  is  also  plotted . 37 

19:  Background  color  represents  the  shear  wave  velocity  (Vs)  from  Moschetti  et  ah, 

(2010)  at  various  crustal  depths,  white  circle  is  the  event  location,  and  black  triangles 
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20:  Waveform  comparison  between  finite-difference  (SW4,  solid  black  line)  and 
frequency  wavenumber  integration  (FK,  dashed  line) . 44 
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21:  The  simulated  3D  data  is  in  black,  the  deviatoric  solution  is  in  green  and  the 
full  solution  is  in  brown.  All  waveforms  are  in  displacement  (cm) . 46 

22:  (a)  Network  Sensitivity  Solutions  color-coded  by  variance  reduction  (VR)  pre¬ 
sented  on  the  Tape  and  Tape  (2012a)  and  Tape  and  Tape  (2012b)  Lune.  The  white 
circle  represents  the  location  of  the  best  full  moment  tensor  solution  in  source-type 
space,  (b)  VR  with  respect  to  source  depth  for  two  different  DC  mechanisms  and 
color-coded  by  source  depth.  The  filled  circles  are  full  moment  tensor  solutions  and 
the  filled  squares  are  deviatoric  moment  tensor  solutions.  The  full  moment  tensor 
solutions  are  also  plotted  in  source-type  space  and  color-coded  by  source  depth  as 


well . 47 

23:  Dashed  lines  are  waveforms  calculated  using  a  ID  Earth  model  and  solid  lines  are 
waveforms  calculated  using  a  3D  Earth  model.  The  waveforms  are  in  displacement 
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24:  (a)  Explosion  moment  tensor  solutions  and  waveform  fits  from  10  to  50  seconds. 
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using  only  waveforms,  only  first  motions  or  combining  waveforms  and  first  motions. 

The  white  circle  represents  the  location  of  the  best  full  moment  tensor  solution  in 
source-type  space . 50 

25:  (a)  Deviatoric  (green)  and  full  (brown)  moment  tensor  solutions  and  waveform 
fits.  The  simulated  data  for  an  explosion  plus  reverse  fault  mechanism  are  in  black. 
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29:  ID  and  3D  synthetic  waveforms  for  a  normal  mechanism  (strike=202; 
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1.  SUMMARY 


In  this  project  (FA9453-13-C-0271)  we  build  on  our  earlier  results  (DE-FC52- 
06NA27324,  Dreger  et  al.,  2008;  Ford  et  al.,  2010;  Ford  et  al.,  2009a;  Ford  et  al.,  2008;  Chiang 
et  al.,  2014,  Nayak  and  Dreger,  2014)  to  investigate  the  effects  of  3D  velocity  structure  on 
seismic  moment  tensor  recover  and  moment  tensor  source-type  discrimination.  In  addition,  we 
have  investigated  the  effect  of  shallow  depth  of  burial  and  the  issue  of  vanishing  traction  at  the 
free-surface  on  fundamental  Green’s  function  amplitudes  and  on  the  recovery  of  the  seismic 
moment  tensor.  Finally  we  have  developed  a  source-type  specific  inversion  that  enables  the 
estimation  of  the  true  maximum  fit  surface  in  the  moment  tensor  source-type  space. 

2.  INTRODUCTION 

The  estimation  of  the  seismic  moment  tensor  using  long-period  regional  distance  three- 
component  waveform  data  has  proved  to  be  an  effective  approach  in  the  discrimination  of 
explosions  from  earthquakes  (Ford  et  al.,  2009a, b;  Chiang  et  al.,  2014).  In  fact,  the  approach  has 
been  found  to  be  able  to  uniquely  discriminate  explosions  from  earthquakes  as  evaluated  in 
source-type  space  (Hudson  et  al.,  1989;  Tape  and  Tape,  2015)  when  long-period  regional 
distance  waveforms  are  combined  with  regional  or  teleseismic  distance  P-wave  first-motion 
observations  (e.g.  Ford  et  al.,  2012;  Chiang  et  al.,  2014).  Sources  of  error  such  as  noise  in  the 
data  and  poor  station  coverage  are  now  routinely  evaluated,  however  the  primary  source  of  error 
is  due  to  the  assumption  of  a  suitable  seismic  velocity  model  from  which  to  calculate  the  needed 
Green’s  functions  for  the  inversion  and  remains  an  important  problem.  Although  it  is  becoming 
common  to  evaluate  uncertainty  due  to  assumed  ID  velocity  models  (e.g.  Ford,  2008;  Chiang  et 
al.,  2014)  the  effect  of  unaccounted  for  3D  velocity  structure  is  an  important  open  problem, 
which  sections  3.3  and  3.4  of  this  report  address.  For  the  western  US  earthquakes  and  the  NNSS 
nuclear  explosions  a  synthetic  investigation  of  the  effects  of  3D  structure  on  the  recovery  of 
seismic  moment  tensors,  particularly  the  isotropic  component  needed  for  source-type 
discrimination  is  presented.  Invoking  source-receiver  reciprocity  (Eisner  and  Clayton,  2001) 
Green’s  functions  for  a  surface  wave  group  and  phase  velocity  constrained  seismic  velocity 
model  (Moschetti  et  al.,  2010)  are  applied  to  NNSS  nuclear  explosions  and  earthquakes  to  assess 
the  recovery  of  the  seismic  moment  tensor  and  its  use  in  source-type  discrimination. 

Also  during  the  period  of  the  project  two  other  important  problems  were  studied.  The 
first,  in  section  3.2,  involves  the  effect  of  the  free-surface  boundary  condition  and  the  vanishing 
nature  of  some  fundamental  moment  tensor  Green’s  functions  for  events  with  shallow  depth  of 
burial.  A  thorough  synthetic  seismogram  analysis  of  this  effect  on  the  recovery  of  seismic 
moment  tensors  and  the  discrimination  potential  was  carried  out.  Shallow  (~1 5m  depth)  mining 
related  explosions  were  then  studied  at  local  distances  to  corroborate  the  synthetic  results,  and  to 
demonstrate  the  ability  to  obtain  moment  tensors  for  these  small  explosions,  discriminate  the 
source-type,  as  well  as  estimate  the  yield  of  the  small  explosions.  The  second,  in  section  3.1, 
involved  the  development  of  a  new  inverse  procedure  for  estimating  seismic  moment  tensors  in 
the  framework  of  the  source-type  projections  (either  Hudson  et  al.,  1989;  Tape  and  Tape,  2015; 
or  others).  This  procedure  called  NSSinv  solves  the  linearized  inversion  for  the  best  fitting 
eigenvectors  for  specific  ratios  of  eigenvalues  (source-type).  The  method  is  fast,  improving 
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substantially  on  the  brute  force  approach  of  Ford  et  al.  (2010),  and  due  to  the  direct  inverse 
nature  of  the  method  it  recovers  the  true  maximum  fit  surface  in  source-type  space.  The  approach 
is  demonstrated  through  inversions  of  events  occurring  in  a  geothennal  environment,  those 
occurring  as  the  result  of  the  collapse  of  a  brine  cavity  in  a  salt  dome,  and  for  a  shallow  mining 
related  explosion. 

3.  TECHNICAL  WORK 
3.1.  Source-Type  Specific  Inversion 

The  mapping  of  the  fit  of  seismic  moment  tensor  (MT)  solutions  in  source-type  space 
helps  to  characterize  uncertainty  and  solution  uniqueness.  Current  practice  relies  on  the  forward 
testing  of  a  distribution  of  randomly  generated  MT  in  source-type  space,  which  is  slow  and  does 
not  necessarily  recover  the  true  maximum  fit  surface.  We  design  an  iterative  damped  least 
squares  inversion  scheme  to  invert  waveforms  and/or  P-wave  first-motions  (FM)  for  best- fitting 
MT  solutions  for  specific  source-types.  An  event  associated  with  the  sinkhole  at  the 
Napoleonville  Salt  Dome,  Louisiana,  an  industrial  quarry  explosion  and  an  earthquake  at  The 
Geysers  geothermal  field,  Northern  California  are  presented  as  examples.  We  find  that  the 
inversion  method  is  more  accurate  and  successful  than  the  random  search  approach  in  recovering 
the  region  of  best- fitting  MT  solutions  or  source-types  and  is  substantially  faster.  The  approach 
also  enables  the  detennination  of  the  best-fitting  MT  for  specified  source-types  such  as  pure 
double-couples,  tensile  cracks  or  explosions,  as  well  as  compound  mechanisms  in  a  single 
numerical  framework.  This  section  (3.1)  has  been  published  in  the  Bulletin  of  the  Seismological 
Society  of  America  (Nayak  and  Dreger,  2015). 

3.1.1  Introduction 

The  2nd  order  general  seismic  moment  tensor  (MT)  is  routinely  used  to  describe  source 
mechanisms  of  seismic  events  (Jost  and  Hernnann,  1989;  Julian  et  al.,  1998;  Aki  and  Richards, 
2002;  Minson  and  Dreger,  2008).  The  general  MT  is  mathematically  defined  to  be  a  symmetric 
3x3  matrix  and  has  six  independent  components  (Mxx,  MYy,  MZ2,  Mxy>  Myz,  Mzz),  which 
describe  the  seismic  source  mechanism  in  terms  of  moments  of  body  force  equivalents,  i.e., 
double  couples  and  linear  vector  dipoles.  Recently,  there  has  been  renewed  interest  in  the 
geometric  representation,  decomposition  and  interpretation  of  the  general  MT  (Vavrycuk,  2011; 
Tape  and  Tape  2012a,  b,  2013;  Vavrycuk,  2015),  which  facilitate  interpretation  of  seismic 
sources.  Source-type  discrimination  and  assessing  confidence  and  uncertainties  in  source-type 
characterization  of  seismic  events  is  of  great  importance,  especially  for  monitoring  of  nuclear 
explosions  (e.g.,  Ford  et  al.,  2009a,  b;  2010,  2012;  Chiang  et  al.,  2014)  and  analysis  of  induced 
seismic  events  (e.g.,  Sileny  et  al,  2009;  Guilhem  et  al.,  2014).  In  this  study,  we  describe  the 
concept  of  a  source-type  specific  MT  inversion.  We  provide  the  expressions  for  a  general  MT  in 
terms  of  its  normalized  eigenvalues,  eigenvectors  and  a  moment  scale  factor.  Then  we  describe 
an  iterative  damped  least  squares  (LS)  inversion  scheme  to  invert  for  best-fitting  eigenvectors 
and  moment  scale  factor  for  an  event  using  its  displacement  waveforms  for  a  given  source-type, 
i.e.,  its  normalized  eigenvalues  that  enables  the  construction  of  the  maximum  goodness  of  fit 
surface  in  the  source-type  space.  We  validate  the  inversion  scheme  by  applying  it  on  synthetic 
and  observed  seismic  waveforms,  and  analyze  the  results.  Finally,  we  also  test  a  method  to 
incorporate  P-wave  first  motion  (FM)  polarities  in  the  inversion. 
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3.1.2  Methodology 


Let  a  given  general  MT,  M,  have  eigenvalues,  Z0  =  \A01  ,AQ2  ,  A03]  and  corresponding 
eigenvectors  e1;  e2  and  e3  (ej  =  [e^  ,  ei2  ,  &i 3]).  The  general  MT  being  real  and  symmetric  has 
eigenvalues  and  eigenvectors  that  are  real  and  orthonormal,  respectively.  Z0  determines  the 
seismic  scalar  moment  and  the  source-type,  and  can  be  nonnalized  to  a  unit  vector,  X,  by  its  L2- 


norm,  m o  (  = 


A0  i 

A  =  JhL  for  i  =  1,2,3  (1) 
mo 

ml  is  a  moment  scale  factor,  defined  as  a  square  so  that  it  remains  non-negative  and 
characterizes  the  absolute  size  of  MT  eigenvalues  (or  of  MT  elements)  independent  of  their  sign. 
Substituting  Aoi  —  At  x  ml  in  equation  (22)  in  Jost  and  Herrmann  (1989),  each  MT  element, 
Mij,  can  be  expressed  as  a  scalar  function  of  the  moment  scale  factor,  normalized  eigenvalues 
and  orthononnal  eigenvectors,  where 

=  ml  (A^e^  +  A2e2ie2j  +  A3e3ie3j )  for  i  =  1,2,3  and  j  =  1,2,3  (2) 

in  which  1,  2  and  3  represent  X,  Y  and  Z  for  the  MT  elements.  The  aim  of  this  study  is  to  assume 
specific  values  of  X,  and  then  solve  for  ml  and  MT  eigenvectors  e1,e2  and  e3,  using  seismic 
waveforms  and/or  polarities.  We  also  discuss  the  utility  of  MT  inversions  whose  solutions  are 
constrained  to  specific  values  of  X. 


Parameterization  of  Eigenvectors 

The  orthononnal  eigenvectors  e1(  e2  and  e3  can  be  expressed  in  terms  of  either  spherical 
( 81,d2,93 )  or  cartesian  (a1,a2,a3,b1,b2)  parameters.  In  this  study  we  develop  the  cartesian 
parameterization  as  we  found  it  to  have  better  convergence  properties.  In  the  electronic 
supplement,  we  provide  the  alternative  parameterization  in  spherical  coordinates  and  show  that  it 
can  be  used  with  the  same  inversion  scheme  to  obtain  nearly  the  same  results,  albeit  with  slower 
convergence  towards  the  best-fitting  solution. 

Assuming  the  five  unconstrained  real  parameters  (a1,a2,a3,b1,b2)  of  the  3D  cartesian 
system,  we  can  define  ex  and  e2 


ei 


[%,  a2,  a3\ 

ri 


(3) 


e2  — 


[a-A  ,  a3b2  ,  -(aA  +  a2b2 )] 

r2 


(4) 


ri 


i=  1 


of 

u-l  , 


r2  =  V (aA)2  +  (a3b2)2  +  (aA  +  a2b2 )2  (5) 
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e3  —  ej  x  e2  (6) 


It  is  important  to  note  that  there  is  innate  numerical  non-uniqueness  in  solutions  of 
orthonormal  eigenvectors  for  a  particular  MT.  There  are  4  combinations  (+e1;  +e2)  of  e1;  e2 
and  e3  that  give  the  same  MT  (equation  2).  In  addition,  the  number  of  independent  parameters 
required  for  defining  3  orthononnal  eigenvectors  is  3.  For  the  cartesian  formulation  in  equations 
(3. 1.3-6),  our  model  is  over-parameterized  in  order  to  allow  the  model  parameters  to  assume  any 
value,  as  long  as  ex  and  e2  are  not  zero  vectors.  As  a  result,  there  can  be  infinite  combinations  of 
(a1,  a2,  a3,  blt  b2)  that  lead  to  the  same  values  of  e1(  e2  and  e3. 

Inverse  Problem  Formulation 


Elements  of  e1(  e2  and  e3  in  equation  (2)  can  be  substituted  by  alr  a2,  a3,  bx  and  b2  to 
obtain  scalar  expressions  of  MT  elements,  Mtj,  in  terms  of  ml,  X  and  the  eigenvector  parameters. 

m' =  f (X,ml,a1,a2,a3,b1,b2)  (7) 

where  m'and  f  are  column  vectors  containing  MT  elements  and  their  scalar  expressions, 
respectively. 


Now  the  MT  elements  in  the  expressions  relating  transverse  (uT),  radial  (uR)  and  vertical 
(uz)  displacements  to  the  MT  and  Green’s  Functions  (GF)  matrix  (G)  in  Minson  and  Dreger 
(2008)  can  be  replaced  by  their  functions  in  terms  of  X,  ml,  a1,  a2,  a3,  blt  and  b2. 

ut,r,z  —  Gm'  (8) 

uT,R,z  =  Gf( X,ml,a1,a2,a3,b1,b2)  (9) 


Equation  (8)  is  valid  for  a  point  source  that  assumes  the  source-time  function  to  be  an 
impulse  or  a  Dirac  delta  function  in  space  and  time  that  is  common  for  all  MT  elements.  For 
constant  and  known  source-type  X ,  we  can  invert  for  the  six  unknown  parameters 
(m0,  a1,  a2,  a3,  blt  b2).  Since  we  know  the  exact  expressions  for  /  and  its  derivatives,  we  choose 
to  solve  the  problem  using  an  iterative  damped  LS  inversion  scheme  as  follows 


/  df  x^3  df  v-12  df  \ 

du™z  =  G(d^dm° +li^ida‘ +li,1Widb‘)  (10) 


Defining  the  model  parameter  vector  dx 


dx  =  [dmo.da^d^.d^.dbi.db^  (11) 


P  =  G 


■  df 
. dm0 


df 

dax 


df  df  df  df  - 
da2  da3  db1  db2. 


(12) 


duTRZ  =  Pdx  (13) 
dx  =  (PTP  +  /cI)_1PTduTRZ  (14) 

where  k  is  a  damping  parameter.  The  data  goodness  of  fit  parameter  is  the  variance  reduction 
(VR;  expressed  in  %),  a  measure  of  normalized  goodness  of  fit  between  observed  and  synthetic 
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data  (Pasyanos  et  al.,  1996).  The  expressions  of  partial  derivatives  in  equation  (12)  are  provided 
in  Appendix. 

Equation  (9)  is  exact  and  the  user  can  choose  to  solve  it  using  any  appropriate  numerical 
inversion  technique.  This  approach  to  the  problem  also  implies  that  it  would  be  possible  to 
estimate  the  best-fitting  source-type  or  eigenvalues,  for  specific  orientation/eigenvectors,  such  as 
for  a  specific  fault  plane  or  DC  focal  mechanism. 

Inversion  Parameters 

The  iterative  damped  LS  inversion  procedure  depends  on  many  parameters  like  the 
damping  parameter,  initial  model  parameters  and  number  of  iterations.  A  few  trials  (55  15)  with 
different  initial  model  parameters  are  usually  sufficient  for  convergence  to  the  best-fitting  MT 
solution.  We  use  randomly  generated  real  numbers  of  the  order  ~  1  for  initial  values  of 
a1,a2,a3,b1  and  b2.  For  m0,  we  use  initial  values  logarithmically  distributed  over  an  order  of 
magnitude  around  ~  yf M0/1.0el3,  where  M0  is  the  scalar  seismic  moment  (in  N.m)  computed 
from  preliminary  analyses  or  general  MT  inversion  using  the  definition  of  Bowers  and  Hudson 
(1999).  It  is  important  to  note  that  the  cartesian  eigenvector  model  is  over-parameterized  and 
therefore,  only  3  out  of  5  columns  in  P,  i.e.  containing  partial  derivatives  of  f  with  respect  to 
(a1,a2,a3,b1,b2),  are  linearly  independent.  Since  P  is  a  rank  deficient  matrix  (rank(P)  =  4),  a 
non-zero  damping  value  ( k )  must  be  used  for  PrP  to  be  invertible.  We  have  found  from  trial  and 
error  that  k  ~  IE-10  performs  well.  Multiple  solutions  of  (a1,a2,a3,b1,b2)  give  the  same  (e1; 
e2),  and  4  combinations  of  ex  and  e2  (±e1(  +e2)  give  the  same  MT  elements  and  therefore,  the 
same  waveform  fits.  Depending  on  the  initial  values,  the  inversion  proceeds  towards  any  of  these 
solutions,  minimizing  the  least  square  error  between  observed  and  predicted  waveforms.  With 
successive  iterations,  linearly  independent  MT  elements  (which  are  functions  of  model 
parameters  and  fixed  A)  converge  towards  their  best-fitting  values.  We  terminate  the  inversion  if 
VR  changes  less  than  0.01%  over  2  successive  iterations. 

3.1.3  Tests  on  Synthetic  Waveforms 

We  use  synthetic  waveforms  to  evaluate  the  effectiveness  of  the  inversion  procedure.  We 
forward  model  synthetic  three-component  displacement  waveforms  at  seismic  stations  in 
Northern  California  using  equation  (8)  with  100  randomly  generated  MT  assuming  random 
eigenvectors  and  random  values  of  m0  (between  10'"  and  10  ")  for  each  of  4  source-types:,  (1) 
pure  DC  (A  =  [0.7071,0,-0.7071]),  (2)  a  tensile  crack  in  Poisson’s  solid  (v  =  0.25;  X  — 
[0.9045,0.3015,0.3015]),  (3)  explosion  (X  =  [0.5774,0.5774,0.5774]),  and  (4)  pure  CLVD 
(A  =  [0.8165,  —0.4082,  —0.4082]).  For  these  hypothetical  events  and  the  real  seismic  events 
described  in  the  subsequent  sections,  the  GF  were  computed  using  appropriate  velocity  and 
density  models  for  each  region  and  a  frequency-wavenumber  integration  code  based  on  Haskell 
(1964),  and  Wang  and  Herrmann  (1980)  as  provided  in  Herrmann  (2013).  The  frequency- 
wavenumber  integration  method  computes  complete  three-component  seismograms  (including 
near-field  terms)  consisting  of  all  body-wave  and  surface- wave  phases  for  ID  layered  velocity 
models.  Table  1  shows  the  basic  infonnation  on  the  seismic  events  in  this  study.  For  the 
randomly  oriented  hypothetical  seismic  events,  filtered  GF  were  used  to  compute  synthetic 
wavefonns  with  the  randomly  generated  MT.  The  same  GF  were  used  for  inversion  to 
investigate  convergence  properties  of  the  linearized  inversion.  Figure  1  shows  the  recording 
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stations  for  synthetic  test  events  located  at  8  km  depth  and  one  example  source  mechanism.  In 
this  study,  we  decompose  all  MT  into  a  combination  of  isotropic,  DC  and  CLVD  MT  assuming 
same  principal  stress  orientations  for  DC  and  CLVD  MT  (e.g.,  Jost  and  Herrmann,  1989),  and 
compute  their  relative  contribution  to  the  total  M0  (DC,  CLVD  and  ISO  expressed  in  %  in  Figure 
3.1.3b).  For  all  4  source-types,  all  respective  100  iterative  inversions,  MT  were  recovered 
correctly  with  all  MTVR  (model  goodness-of-fit  comparing  the  inverted  MT  solution  with  the 
actual  MT  used  to  compute  synthetic  data)  and  wavefonn  VR  greater  than  99.5%  and  99.7%, 
respectively. 


(a)  Northern  California  (b) 

JRSC 

PKD 

BKS 

CMB 

KCC 

ORV 

124*W  122°W  120'W 


langentiai  Haaiai  a  ver 

AF~ /yW yJY 

R  =  99.6  km  Az  =  31 6=  Dmax  »  1  4e-04  cm  a 

'y\A'~ —  /\A— —  Vv 


1 37°  Dmax  =  1  3e-04  cm 


AA/V - '\l\p — 


Vertical  Depth  =  8.0  km 

_ Strike  =  96°  ;  191° 

t _ ,  Rake  =  169°  ;  29° 

30  s  Dip  =  61° ;  81° 

M0  =  3.35e+15  N.rr 

- M^  =  4.29 

DC  =  100% 

CLVD  =  0% 

ISO  =  0% 


332°  Dmax  =  9.6e-05  cm 


34°  Dmax  =  1  8e-04  cm 


=  71°  Dmax  =  1 ,5e-04  cm 


=  359°  Dmax  =  1  8e-04  cm 


jrI?s 

•  PKD 


Figure  1.  (a)  Map  shows  the  epicenter  (black  star)  of  a  hypothetical  seismic  event,  whose 
synthetic  waveforms  at  Berkeley  Digital  Seismic  Network  stations  (black  triangles)  were  used 
for  evaluation  of  the  inversion  procedure  proposed  in  this  study,  assuming  different  source-types, 
(b)  Three-component  synthetic  displacement  waveforms  (0.02-0.05  Hz)  for  one  example  random 
MT  solution  assuming  a  DC  source-type.  R  =  epicentral  distance,  Az  =  azimuth,  and  Dmax  = 
maximum  displacement  amplitude  at  a  station.  Beach-ball  represents  lower  hemisphere  P-wave 
radiation  pattern. 
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Table  1.  Event  infonnation  and  parameters  of  MT  inversion  for  the  three  events  in  this  study. 


Hypothetical 

Events 

Event  TE1,  NSD 
Sinkhole 

Event  TE2, 
HUMMING 
ALBATROSS 

Event  TE3,  The 
Geysers,  Northern 
California 

Reference 

- 

Nayak  and  Dreger 
(2014,  2015) 

Chiang  et  al.  (2016) 

Boyd  et  al.,  (2015) 

Date  (yymmdd) 

- 

20120801 

- 

20110301 

Origin  time 
(hh:mm:ss.ss) 

- 

20:52:38.50 

- 

02:19:47.01 

Hypocenter- 
Longitude  (°E), 
Latitude  (°N), 
Depth  (km) 

-121.464, 
36.755,  8.0 

-91.1422,  30.0112, 
0.47 

Depth  =  9  m 

-122.8200, 

38.8153,3.5 

Mw 

- 

1.36 

1.89 

4.50 

Type  of 

waveforms  used 

Displacement 

Displacement 

Velocity 

Displacement 

Filter 

(Butterworth) 

0.02-0.05  Hz, 
p  2  n  2 

0.1-0. 2,  or  0.1-0. 3 
Hz,  p  1  n  4 

1. 2-2.0  Hz  p  2  n  2 

0.02-0.05  Hz,  p  2 
n  4 

Recording  network 

Berkeley 

Digital 

Seismic 

Network 

US  Geological 
Survey  Temporary 
Network 

T  emporary 
broadband  and  short 
period 

seismometers 

Berkeley  Digital 
Seismic  Network, 
Northern 

California  Seismic 
Network, 
Lawrence 
Berkeley  National 
Laboratory  Short 
Period  Network  at 
The  Geysers 

Velocity  Model 

gil7  (Pasyanos 
et  al.,  1996) 

Nayak  and  Dreger 
(2014) 

Saikia  et  al.  (1990) 

gil7,  SoCal 
(Dreger  and 
Helmberger,  1993) 

Number  of  stations 
used  for  waveform 
inversion 

6 

5 

5 

11 

Distance  range  of 
stations  for 
waveform 
inversion  (km) 

100-311 

0.4-1. 4 

1.2-4. 3 

61-230 

Number  of  first 
motion  P-wave 
polarities 

- 

6 

16 

173 

Inverse  distance 
weights  for 
waveform-only 
MT  Inversion 

yes 

no;  inverse  variance 
weights  used 
instead 

yes 

yes 
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3.1.4  Tests  on  a  Real  Event 


Numerous  seismic  events  were  associated  with  the  development  of  a  sinkhole  at  the 
western  edge  of  the  Napoleonville  salt  dome  (NSD),  Louisiana  in  2012.  MT  inversion  of  these 
events  using  a  grid-search  approach  and  separate  preliminary  ID  velocity  models  for  the  salt 
dome  and  the  surrounding  sediment  sequence  yielded  large  isotropic  volume-increase 
components  in  the  MT  solutions  (Nayak  and  Dreger,  2014).  The  centroid  locations  of  these 
events  were  found  to  be  at  the  western  edge  of  the  salt  dome  at  ~  470  m  depth.  Here  we  use  the 
iterative  damped  LS  inversion  scheme  to  estimate  the  MT  solution  of  one  of  these  events  (event 
TE1  in  Nayak  and  Dreger  [2014])  assuming  various  source-types.  In  this  study,  we  consider  pure 
opening  cracks  assuming  sources  in  soft  sediments  (v  ~  0.43)  and  salt  (v  ~  0.24),  a  pure  DC 
source  and  a  pure  explosion.  A  corresponding  to  these  four  sources  are  indicated  in  Figure  2.  For 
the  MT  inversion  of  this  event  and  other  real  seismic  events  described  subsequently,  three- 
component  broadband  velocity  waveforms  were  first  corrected  for  instrument  response, 
integrated  to  displacement  and  then  filtered  in  a  low-frequency  pass  band  appropriate  for  each 
seismic  event.  The  horizontal,  east-west  and  north-south  components  were  rotated  to  radial  and 
transverse  components.  GF  were  filtered  in  the  same  pass  band  as  the  observed  waveforms. 


(a)  Tangential  Radial  Vertical 
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R  =  0  4  km  Az  =  243°  Drnax't  1 ,5e-04  cm  VR  =  82% 
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R  =  0.4  km  Az  =  96"  of  ax  =  7.2e-05  cm  VR  =  34% 

A  =  [0.6845,  0.5155,  0.5155] 

[  a,  a ,  a3  b,  b 3  m*]  =  [  2.13,  -3.83,  -0.77,  -0.79, 1.10,  0.0196  ] 
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A  =  [  0.9150,  0.2853,  0.2853  ] 

[  a,  a2  a3  b,  b2  m/1  =  (  5.27,  -8.42,  -4.87,  -1.10,  0.37,  0.0060  ] 


(C)  Tangential 
LA01 


M0  =  6.74e+10  N.m 
M^*  1.16 
DC  =  0% 

_  CLVD  =  0% 

6.25  s  ISO  =  100% 

VR  =  32.3% 


(d) 
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-A/1  Aac/v 


VPf  =  43% 


LA08 


A  =  [  0.5774,  0.5774,  0.5774  ] 
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Figure  2.  Best- fitting  MT  solutions  of  event  TE1  described  by  normalized  eigenvalues  (A), 
eigenvector  parameters  (a1,a2,a3,b1,b2)  and  seismic  moment  scale  factor  (m^)  for  various 
source-types:  (a)  tensile  crack  in  sediment,  (b)  tensile  crack  in  salt,  (c)  pure  explosion  and  (d)  a 
DC.  The  solid  and  dashed  lines  are  observed  and  synthetic  displacement  waveforms, 
respectively.  Meaning  of  other  symbols  is  same  as  in  Figure  lb.  Station-specific  R,  Az  and 
Dmax  are  same  in  all  subplots.  Waveforms  were  filtered  in  the  pass  band  0. 1-0.3  Hz  for  stations 
LA01,  LA02,  LA03  and  LA09,  and  in  0. 1-0.2  Hz  for  station  LA08.  Final  values  of  (aj,  a2,  a3,  b]% 
b2)  are  non-unique  (see  Figure  A2,  Appendix). 
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The  best  fitting  MT  solutions  and  waveform  fits  for  various  source-types  for  event  TE1 
are  shown  in  Figure  2.  A  DC  source,  an  explosion  and  a  pure  crack  in  salt  fit  the  waveforms 
poorly  at  VR  31.5%,  32.3%  and  45.7%,  respectively.  However,  a  pure  crack  in  sediments  fits  the 
waveforms  well  at  66.6%,  which  is  slightly  lower  than  VR  for  the  full  MT  solution,  68.3% 
(Nayak  and  Dreger,  2014).  The  uniform-phase  nature  of  long-period  S-waves  observed  in  radial 
components  at  all  stations  favors  a  volumetric  source,  whereas  the  strong  variation  of  amplitudes 
with  azimuth,  and  the  presence  of  SH  waves  favors  a  tensile-crack-type  source  rather  than  a 
spherically  symmetric  explosion.  The  strike  of  the  tensile-crack  plane  for  a  crack  in  sediments 
(30°  or  208°)  agrees  very  well  with  the  strike  of  the  DC  component  in  the  full  MT  solution  (17° 
or  232°)  and  with  the  strike  of  a  shear-tensile  source  (crack  +  DC)  in  sediments  (23°  or  225°) 
estimated  in  Nayak  and  Dreger  (2014).  To  verify  our  results,  we  estimated  best  fitting  MT 
solutions  for  these  source-types  using  a  grid  search,  and  found  the  results  agreed  very  well  with 
our  inversion  results.  With  this  method  separate  constrained  and  linearized  LS  inversion 
formulations  or  grid  search  formulations  for  common  source-types  like  explosions  (Ford  et  al., 
2009b),  DC  (Herrmann  et  al.,  2011),  crack  and  pipe  (Nakano  and  Kumagai,  2005;  Minson  et  al, 
2007)  sources  can  be  replaced  by  a  single  mathematical  fonnulation  and  inversion  procedure 
where  the  user  is  required  to  only  specify  appropriate  nonnalized  eigenvalues  for  a  particular 
source-type.  Moreover,  the  grid-search  approach  to  estimate  the  best  fitting  values  of  M  0  and 
fault  plane  parameters  (e.g.,  <p  [strike],  (  [rake]  and  8  [dip]  for  pure  DC  MT  solutions)  is 
generally  time  consuming  (Herrmann  et  al.,  2011). 


3.1.5  Application  to  Maximum  Fit  Surfaces  in  Source-Type  Space 


The  maximum  fit  surface  in  source-type  space  has  been  called  the  Network  Sensitivity 
Solution  (NSS)  in  Ford  et  al.  (2010)  due  to  its  ability  to  assess  recovery  of  source-type 
infonnation  under  changing  network  topology.  The  NSS  is  also  used  to  assess  confidence  in 
source  mechanisms  of  seismic  events  that  are  obtained  from  MT  inversion  with  respect  to  its 
source-type  (or  nonnalized  eigenvalues).  The  NSS  compares  fits  between  observed  and  synthetic 
waveforms  for  a  large  population  (usually  on  the  order  of  tens  of  millions)  of  MT  covering  the 
entire  source-type  space  (Ford  et  al.,  2010).  MT  are  assembled  from  random  populations  of 
eigenvalues  and  orthonormal  eigenvectors.  The  eigenvalues  are  randomly  drawn  from  a 

1 1  UT,R,Z  1 1 


population  of  real  numbers  unifonnly  distributed  between  —nL  and  nL  where  L  — 


IIGII 


is  a 


factor  representing  absolute  size  of  the  eigenvalues  and  the  value  of  n  (usually  >  7)  specifies  the 
range  of  the  eigenvalues  with  respect  to  L.  Since  the  NSS  depends  on  uT  R  Z  and  G,  it  takes  into 
account  the  station  distribution,  frequency  content  of  waveforms,  and  data  quality  and  quantity 
for  a  given  MT  inversion  scenario.  It  generates  a  distribution  of  VR  in  source-type  space  that 
allows  us  to  identify  the  uniqueness  of  the  source-type  obtained  from  MT  inversion,  and  the 
existence  of  possible  tradeoffs  commonly  observed  in  nuclear  explosions  (Ford  et  al.,  2010; 
Chiang  et  al.,  2014).  The  eigenvalues  are  used  to  compute  the  two  source-type  parameters:  e  and 
k  (Hudson  et  al.,  1989).  In  the  Hudson  plot  (Hudson  et  al.,  1989)  the  horizontal  axis  plots  the 
ratio  of  the  deviatoric  eigenvalues  (c)  and  the  vertical  axis  plots  the  relative  isotropic  component 
(k).  For  each  coordinate  (c,  k),  the  best  VR  value  from  all  MT  corresponding  to  a  small  area 
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around  that  coordinate  is  selected  and  plotted  to  generate  the  NSS  (Ford  et  al.,  2010). 


x  -CLVD 
[1,1. -2] 


-Dipole 
[0.0, -1] 


-Crack  [-1  ,-1  ,-3] 


Implosion  [-1,-1  ,-1 


Figure  3.  Grid  of  7457  unique  normalized  eigenvalues  (black  ‘+’  signs)  or  unique  source-types 
on  (a)  the  fundamental  Lune  (Tape  and  Tape,  20 12a, b),  and  (b)  on  the  Hudson  source-type  plot 
(Hudson  et  al.,  1989).  Black  crosses  are  positions  of  major  theoretical  source-types  shown  with 
their  un-normalized  eigenvalues. 


The  forward  modeling  of  synthetic  waveforms  for  millions  of  MT  required  to  produce  the 
NSS  is  computationally  intensive  and  time  consuming.  Moreover,  considering  a  6D  parameter 
space,  the  population  of  random  MT  needs  to  be  very  large  in  order  to  get  some  of  the  MT  close 
to  the  true  best-fitting  solutions  for  a  particular  source-type,  which  can,  at  best,  be  approximate. 
However,  we  can  use  our  source-type  specific  MT  inversion  method  to  compute  the  NSS  with 
greater  accuracy  and  efficiency.  We  first  construct  a  grid  on  the  fundamental  Lune  of  the 
normalized  eigenvalue  sphere  using  equation  (20)  of  Tape  and  Tape  (2012a),  which  is 
reproduced  in  equation  (15)  with  the  radial  coordinate  set  to  unity. 


\h] 
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We  keep  the  grid  spacing  in  co-latitude  (1  fixed  at  1.8°  and  decrease  the  grid  spacing  in 
longitude  y  linearly  from  ~  1°  at  the  poles  to  ~  0.6°  at  latitudes  ±  65°,  keeping  it  fixed  at  ~  0.6° 
around  the  equatorial  region  (Figure  3a).  This  grid,  comprising  of  7457  unique  coordinates, 
representing  unique  sets  of  normalized  eigenvalues  or  unique  source-types  (A),  is  projected  on 
the  Hudson  plot  (Figure  3b).  The  eigenvalue  sets  are  arranged  in  a  sequence  having  continuity  in 
source-type  space.  For  each  X,  we  apply  our  iterative  damped  LS  inversion  scheme  to  estimate 
best-fitting  values  of  m0,a1,a2,a3,b1  and  b2.  To  estimate  initial  model  parameter  values  for 
each  X,  we  construct  an  initial  population  of  ~50  solutions  of  m0,  av  a2,  a3,  b1  and  b2  using 
values  from  inversion  results  of  the  previous  X  (thus  exploiting  continuity  of  MT  solutions  in 
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source-type  space)  and  additional  values  of  m0  logarithmically  distributed  over  two  orders  of 
magnitude.  This  set  of  solutions  is  used  to  forward  model  synthetic  waveforms  and  the  solution 
that  returns  the  best  VR,  is  used  as  the  initial  model.  Sequentially  repeating  this  process  for  all  X 
in  the  grid  generates  best-fitting  eigenvectors  and  moment  scale  factors  for  all  X,  yielding  the 
best-fitting  VR  surface  covering  the  fundamental  Lune  or  the  Hudson  space.  We  also  skip  over 
one  or  two  eigenvalue  sets  in  the  grid  if  VR  at  the  previous  X  is  between  -20%  and  30%,  or  less 
than  -20%,  respectively.  This  increases  the  overall  speed  at  the  cost  of  VR  resolution  in  those 
regions  of  the  source-type  grid  that  fit  the  waveforms  poorly. 

We  compute  NSS  for:  (1)  event  TE1  of  the  NSD  sinkhole  sequence  (Nayak  and  Dreger, 
2014)  described  above,  (2)  event  TE2,  a  chemical  explosion  for  industrial  applications  (Chiang, 
A.,  et  al,  2015),  and  (3)  event  TE3,  an  earthquake  at  The  Geysers  geothermal  field,  Northern 
California  (Table  1). 

TE2  is  shot  800  in  HUMMING  ALBATROSS,  an  industrial  quarry  blast  experiment 
involving  a  set  of  chemical  explosions  at  very  shallow  depths  (-10  to  15  m)  and  recorded  at 
distances  up  to  several  kilometers  away  (Chiang  et  al,  2016).  Data  of  these  explosions  have  been 
analyzed  in  detail  to  study  MT  solutions  in  different  frequency  bands,  depth  dependence  of  MT 
solutions  and  source-type  discrimination  using  both  wavefonns  and  FM  polarities  (Chiang  et  al, 
2016). 


TE3  is  an  earthquake  at  The  Geysers  geothennal  field  in  Northern  California  that  has 
witnessed  a  considerable  increase  in  number  of  small  magnitude  earthquakes  (Mw  1.5  to  4.0)  in 
response  to  steam  production  and  water  injection  for  reservoir  recharge  since  the  1960s. 
Anomalous  isotropic  components  have  been  detected  in  full  MT  solutions  of  many  Mw  >3.5 
earthquakes  at  The  Geysers.  The  event  we  study  (TE3)  is  a  Mw  4.5  earthquake  on  1  March  2011, 
that  has  a  well  constrained  30%  isotropic  volume-increase  component  in  its  full  MT  solution 
(Boyd  et  al,  2015). 

The  NSS  results  for  the  three  events  are  shown  in  Figure  4.  MT  solutions  that  produce  the 
best  fits  for  the  sinkhole  event  TE1  are  tightly  clustered  in  a  region  between  theoretical 
explosions  and  tensile  cracks,  quite  far  away  from  theoretical  deviatoric  mechanisms  and  the 
expected  closing  cracks  for  a  collapse  process,  which  produce  fits  only  up  to  VR  <  35%  (Figure 
4a).  Comparison  of  the  shape  of  the  NSS  in  Figure  4a  computed  using  0. 1-0.3  Hz  waveforms, 
with  NSS  in  Figure  14a  in  Nayak  and  Dreger  (2014)  computed  using  0. 1-0.2  Hz  waveforms, 
demonstrates  that  the  constraints  on  source-type  are  stronger  for  NSS  computed  using 
waveforms  in  a  broader  frequency  pass  band. 
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Figure  4.  NSS  of  the  three  events  using  low-frequency  displacement  waveforms:  (a)  TE1,  (b) 
TE2  and  (c)  TE3.  Left  panels  (NSS  Inversion)  show  NSS  computed  using  the  inversion  approach 
in  this  study.  Right  panels  (NSS  Random)  show  NSS  computed  using  randomly  generated  80 
million  MT.  The  contours  and  shading  represent  absolute  values  of  VR  (%)  whereas  NSS  plots 
in  other  studies  usually  show  nonnalized  VR  (e.g.,  Guilhem  et  al.,  2014;  Nayak  and  Dreger, 
2014;  Chiang  et  al,  2014;  Boyd  et  al.,  2015).  For  each  event,  the  VR  scale  is  same  for  both  NSS 
plots  (left  and  right)  to  enable  better  comparison.  Black  crosses  are  positions  of  major 
theoretical  source-types.  For  each  event  the  white  star  is  the  position  of  the  best-fitting  full  MT 
solution  from  a  time-domain  full  MT  inversion  of  waveforms.  In  each  plot  the  white  circle  is  the 
source-type  corresponding  to  the  maximum  VR  recovered  by  each  NSS  (VRMax  in  the  lower 
left  comer).  The  color  version  of  this  figure  is  available  only  in  the  electronic  edition. 
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For  event  TE2,  which  is  a  chemical  explosion,  we  observe  that  both  DC  and  crack-like 
volume-increase  MT  solutions  produce  similar  quality  of  waveform  fits  at  VR  82%-85%  (Figure 
4b),  which  differs  from  common  behavior  of  the  explosion  NSS.  Usually  CLVD  sources  with  the 
vertical  axis  in  compression  provide  similar  quality  of  fits  to  explosion  waveforms  because  an 
isotropic  volume-increase  MT  solution  and  the  vertically-oriented  negative  CLVD  both  have  an 
isotropic  Rayleigh  wave  excitation  and  no  Love  wave  excitation  (Ford  et  al.,  2009b,  2010, 
2012).  The  very  shallow  depth  of  event  TE2  (~  9  m)  causes  the  vertical  dip-slip  (DS) 
components  of  the  GF,  i.e.,  TDS,  RDS,  ZDS  (Minson  and  Dreger,  2008),  to  be  vanishingly  small 
due  to  free-surface  vanishing  traction.  This  coupled  with  the  fact  that  the  source  process  also 
excited  significant  SH  waves  introduces  spurious  M  and  M  components  in  the  MT 
solution  leading  to  the  vertical  dip-slip  nature  of  best-fitting  solutions  near  the  DC  region 
(Chiang  et  al.,  2016).  This  is  further  confirmed  by  inverting  wavefonns  of  event  TE2  assuming  a 
pure  DC  source,  which  returns  a  vertical  dip-slip  MT  solution  with  Mw  1.41,  VR  82.5  %  and  one 
of  the  fault  plane  solutions  defined  as tp  ,  (  and  8  =  345°,  83°  and  84°,  respectively. 


The  best- fitting  region  of  MT  solutions  for  event  TE3  (VR  >  79%)  covers  a  large  area  on 
the  Hudson  plot  (Figure  4c),  and  therefore  its  source-type  is  poorly  constrained  by  the  NSS 
computed  using  low-frequency  displacement  wavefonns  alone.  Its  shape  is  different  from  that  of 
a  typical  earthquake  NSS  (Ford  et  al.,  2010;  Chiang  et  al.,  2014),  which  can  be  used  to  flag 
unusual  events  that  warrant  further  investigation.  Like  many  other  events  at  The  Geysers 
geothennal  field  the  best-fitting  full  MT  solutions  of  event  TE3  show  primarily  positive  isotropic 
components  (Boyd  et  al.,  2015). 


For  comparison,  we  also  compute  the  NSS  for  events  TE1,  TE2  and  TE3  using  the 
forward-modeling  approach  (e.g.  Ford  et  al.,  2010)  with  a  population  of  80  million  MT  (i  =  8) 
for  each  event  (right  panels  in  Figure  4).  For  the  three  events,  we  compare  1)  the  VR  for  the  best- 
fitting  randomly  generated  MT  solution  (VRmax  in  NSS  Random  plots  in  Figure  4),  2)  the  VR 
for  the  best-fitting  MT  solution  obtained  using  our  iterative  damped  LS  inversion  method 
(VRmax  in  NSS  Inversion  plots  in  Figure  4),  and  3)  the  VR  of  the  best-fitting  full  MT  solution 
computed  from  time-domain  full  MT  inversion  of  wavefonns  (Minson  and  Dreger,  2008).  For 
the  three  events  the  best- fitting  random  NSS  VR  are  65.5%,  83.3%  and  80.2%,  respectively.  For 
the  NSS-inversion  method,  they  are  68.2%,  85.4%  and  80.2%,  respectively.  Best-fitting  full  MT 
solution  (Nayak  and  Dreger,  2014;  Chiang  et  al.,  2016;  Boyd  et  al.,  2015)  VR  are  68.3%,  85.4% 
and  80.2%,  respectively.  The  white  circles  and  stars  in  Figure  4  show  the  respective  best-fit 
solutions. 


Even  with  80  million  MT,  the  best-fitting  randomly  generated  MT  solutions  for  events 
TE 1  and  TE2  fit  their  respective  waveforms  at  VR  lower  than  the  VR  of  the  best-fitting  full  MT 
solutions.  The  best  values  of  VR  recovered  by  the  NSS  computed  using  our  iterative  damped 
LS  inversion  approach  are  close  to  or  the  same  VR  as  the  best-fitting  full  MT  solutions  for  all 
three  events.  We  are  also  able  to  recover  the  source-type  of  the  best-fitting  full  MT  solution  as 
seen  in  the  overlap  of  the  white  star  and  the  white  circle  on  the  NSS-inversion  plots  for  all 
three  events.  Overall,  we  observe  that  our  inversion-based  MT  solutions  produce  better 
wavefonn  fits  than  MT  solutions  from  a  population  of  randomly  generated  solutions  by 
VR  ~  0-3%  over  a 
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substantial  area  on  the  Hudson  plot  for  all  three  events. 

Computing  the  NSS  by  estimating  the  best-fitting  MT  solution  at  each  grid  point  on  the 
Hudson  plot  also  gives  us  the  flexibility  to  make  the  grid  coarser  or  liner  depending  on  the 
purpose  of  our  analysis  since  the  number  of  grid  points  affects  the  computation  time.  We  can 
also  choose  to  evaluate  only  a  portion  of  the  source-type  space.  For  example,  we  can  compute 
the  NSS  only  for  source-types  with  positive  sum  of  eigenvalues  for  the  purpose  of  explosion 
monitoring.  For  event  TE2,  which  has  a  total  number  of  7500  wavefonn  samples  (5  stations  x  3 
components  x  500  samples  per  time  series),  a  Fortran  90  code  supported  by  LAPACK  and  Xcode 
running  in  Mac  OS  X  10.8.5  on  a  3rd  generation  Intel  2.6  GHz  i7-3720QM  processor  took  little 
over  3  minutes  in  computing  the  NSS  using  the  inversion-based  approach.  This  opens  up  the 
possibility  of  near-real  time  source-type  confidence  analysis  of  seismic  events  after  a  general  MT 
solution  has  been  computed  (and  possibly  reviewed  by  an  analyst),  using  minimal  computational 
resources.  In  comparison,  the  same  computer  took  ~60  minutes  for  forward  modeling  wavefonns 
with  80  million  random  MT. 

The  method  to  detennine  the  NSS  from  a  random  population  of  eigenvectors  and 
eigenvalues  requires  using  a  uniform  distribution.  Vavrycuk  (2015)  has  shown  that  care  is 
needed  in  the  use  of  various  MT  norms  applied  to  solution  distributions  in  different  source-type 
projections  because  a  uniform  distribution  in  one  projection  may  not  be  uniform  in  a  different 
projection.  Moreover,  distributions  and  uncertainties  in  MT  elements  are  projected  on  source- 
type  plots  in  complicated  ways  on  different  sections  of  the  source-type  plots  (Vavrycuk,  2015). 
However,  the  NSS-inversion  approach  does  not  depend  on  any  assumption  on  the  distribution  of 
MT  eigenvalues.  Instead  the  NSS-inversion  systematically  finds  the  best-fitting  MT  solution 
corresponding  to  each  source-type  (normalized  eigenvalues)  on  a  predefined  source-type  grid 
sufficiently  covering  the  entire  source-type  space  or  a  specific  section  under  investigation,  for 
any  source-type  projection. 

3,1.6  Incorporating  First  Motions  in  the  NSS 

To  better  constrain  the  source  radiation  patterns  of  shallow  explosion  events,  Ford  et  al. 
(2012)  introduced  the  inclusion  of  teleseismic  P-wave  FM  polarities  in  the  NSS  approach  for 
better  sampling  of  the  entire  focal  sphere.  Since  waveforms  used  in  MT  inversion  are  usually  of 
a  low  frequency  nature,  inclusion  of  P-wave  FM  polarities  from  regional  or  teleseismic  distance 
stations  provides  independent  high  frequency  information.  Inversion  of  MT  solutions  from  P- 
wave  FM  polarities  using  the  derivative-based  scheme  that  we  have  implemented  for  wavefonns 
is  problematic  as  the  sign  function  is  not  differentiable  at  one  point,  and  its  derivatives  are  zero 
everywhere  else  in  the  domain.  Therefore,  we  employ  a  continuous  function  to  approximate  the 
sign  function. 

P-wave  FM  polarity  data,  uFM,  which  is  -1  (down  or  dilation),  0  (zero)  or  +1  (up  or 
compression),  can  be  written  as 

1^0  z 

ufm  —  1 - r  for  uoz  ^  0  and  uFM  —  0  otherwise  (16) 

I  uoz  I 

where  uozis  first  arriving  P-wave  vertical  component  displacement  amplitude  that  depends  on 
MT  elements,  as  well  as  azimuths  and  takeoff  angles  to  individual  stations  (expression  for  far- 
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field  P-wave  radiation  in  equation  4.29  in  Aki  and  Richards,  2002).  uFM  can  be  approximated  as 

Ur\  Y 

uFM=  .  (17) 

Ve  +  uoz2 

where  e  is  a  small  positive  number.  Figure  5  compares  the  sign  function  and  its  approximation 
with  e  —  lxl O'4  used  in  this  study,  which  makes  it  continuous  and  differentiable.  uoz  can  be 
expressed  as 

uoz  =  G'm'  =  G'  f  (ml,  X,  ev  e2,  e3)  =  G'  g(X,  av  a2,  a3,  blt  b2)  (18) 

where  G'  is  a  function  of  takeoff  angles  and  azimuths.  Since  we  are  considering  only  the 
polarities  of  the  displacement,  we  can  drop  the  moment  scale  factor,  ml,  as  a  variable  (or 
g  =  f /ml).  Therefore,  for  constant  and  known  source-type  X,  FM  polarity  information  can  be 
expressed  as  an  approximate  scalar  function  of  five  independent  parameters.  Similar  to  the 
procedure  applied  to  waveforms,  we  can  set  up  an  iterative  damped  LS  inversion  scheme  to 
estimate  best-fitting  values  of  (a1;  a2,  a3,  bt,  b2)  for  a  given  X  with  duT  R  Z  and  P  in  equations 

(10-14)  replaced  by  duFM  and  ^-1,  respectively,  with  a  different 

v  5u0z  L  9fli  9^2  9a3  9 b1  9 b2\ 

damping  parameter  k'. 


x 


Figure  5.  Figure  showing  the  true  sign  function  sign(x )  (equation  16),  its  approximation  used  in 
this  study,  uFM(x)  (equation  17),  and  its  normalized  derivative  du™(x\  where  x  —  uoz. 


The  derivatives,  are  specific  to  each  polarity  data  point.  From  Figure  5,  it  is  obvious 

duoz 

that  the  derivatives  will  be  non-zero  only  for  data  points  close  to  the  P-wave  nodes.  Therefore, 
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for  the  FM  polarity  based  inversion  to  proceed,  there  must  be  a  sufficient  amount  of  data,  such 
that  some  data  points  are  always  close  to  zero  and  their  derivatives  are  non-zero.  So,  this 
derivative-based  scheme  for  inversion  of  MT  solutions  works  best  on  large  FM  polarity  data  sets. 
The  quantity  of  FM  polarity  data  required  will  depend  on  the  takeoff  angles,  azimuths  and  the 
source  mechanism  itself.  We  fixed  the  value  of  e  at  IE-4  from  trial  and  error.  Greater  values  of  e 
make  the  inversion  more  stable,  but  the  derivatives  are  more  approximate,  making  the  final 
solution  deviate  from  the  true  solution. 

We  use  this  inversion  scheme  to  construct  a  first-motion  polarity  based  NSS  similar  to 
that  implemented  for  waveforms.  We  apply  this  scheme  on  FM  polarities  of  event  TE3  that  has 
the  largest  amount  of  first-motion  data  among  the  three  events  (173  FM  polarities  collected  from 
the  stations  of  three  networks  in  Northern  California).  For  all  events  in  this  study,  FM  polarities 
were  picked  by  analysts  and  assigned  weights  on  a  scale  of  1  to  3  based  on  their  quality.  We  also 
increase  number  of  initial  models  for  each  A-specific  inversion  to  100  and  we  proceed  with  the 
inversion  with  30  best  initial  models  after  comparing  the  initial  fits.  From  trial  and  error,  k'  = 
IE-2. 


Prior  to  inversion  of  event  TE3  polarity  data,  the  inversion  scheme  was  applied  to 
synthetic  P-wave  FM  polarity  data  sets  for  the  100  random  MT  of  pure  DC  and  pure  CLVD  that 
were  used  in  the  synthetic  wavefonn  tests.  The  azimuths  and  takeoff  angles  in  event  TE3  dataset 
were  used  to  generate  synthetic  P-wave  FM  polarities.  For  both  source-types,  the  inverted  MT 
solutions  are  able  to  fit  the  polarity  data  perfectly  (VR  =  100%)  for  almost  all  (98  out  of  100) 
random  MT,  while  the  final  data  VR  for  the  remaining  MT  is  >  98%.  As  there  can  be  multiple 
similar  eigenvectors  or  orientations  that  fit  a  polarity  dataset  equally,  the  MTVR  corrected  for 
ml  varies  from  -91%  to  99.9%. 

The  results  for  event  TE3  are  shown  in  Figure  6.  While  the  FM  NSS  constrains  the 
mostly  likely  source-type  to  be  close  to  DC  or  deviatoric,  the  maximum  VR  being  only  51% 
indicating  that  there  are  no  source-types  that  fit  all  of  event  TE3  FM  data  well.  This  might  be 
indicative  of  a  complicated  initial  rupture  process,  errors  in  FM  picks,  or  possible  errors  in  depth 
or  takeoff  angles.  The  misfit  between  observed  and  predicted  polarities  (beach-ball  in  Figure  6) 
is  largely  contributed  by:  (1)  the  anomalous  polarities  of  Pn  waves  (polarities  with  steep  equal 
takeoff  angles  in  northeast  and  southeast  quadrants)  that  are  usually  associated  with  uncertainties 
owing  to  their  emergent  nature;  (2)  and  the  section  of  the  focal  sphere  towards  south-southeast 
with  both  positive  and  negative  polarities.  Figure  6  also  shows  the  NSS  computed  from  forward 
modeling  FM  polarities  using  a  population  of  80  million  randomly  generated  MT.  Comparison  of 
the  two  NSS  plots  in  Figure  6  shows  that  our  approximate  inversion  method  works  well  on  large 
FM  polarity  datasets.  It  can  possibly  be  used  along  with  FM  polarity  analysis  software  like 
FPFIT  (Reasenberg  and  Oppenheimer,  1985)  and  HASH  (Hardebeck  and  Shearer,  2002)  to 
identify  anomalous  events  in  near-real  time.  The  approximation  to  the  sign  function  (equation 
17)  doesn’t  affect  the  values  of  VR  as  the  approximation  is  used  to  only  compute  the  derivatives 
whereas  the  actual  sign  function  is  used  to  forward  model  synthetic  FM  polarities  to  compute 
VR. 
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Figure  6.  NSS  of  event  TE3  using  P-wave  first  motion  polarities.  ‘+Dipole’  and  ‘+Crack’  are 
abbreviated  to  ‘+D’  and  ‘+C’,  respectively  in  the  left  panel.  Beach-balls  represent  P-wave  lower 
hemisphere  radiation  pattern  predicted  by  MT  solution  corresponding  to  the  maximum  VR 
recovered  by  each  NSS  (white  circle).  In  the  beach-balls,  black  crosses  and  circles  represent 
observed  positive  and  negative  P-wave  FM  polarities,  respectively;  the  size  of  the  polarity 
symbols  is  scaled  by  their  quality  weight  (1,  2  or  3).  ‘P’  and  ‘T’  indicate  pressure  and  Tension 
axes,  respectively.  Explanation  of  other  symbols,  shading  and  contours  in  NSS  plots  is  same  as 
in  Figure  4.  The  color  version  of  this  figure  is  available  only  in  the  electronic  edition. 


In  order  to  test  if  the  inversion-based  NSS  approach  works  for  a  smaller  quantity  of  FM 
polarity  data,  we  randomly  selected  two  sub-populations  containing  only  30  and  7  FM  polarities, 
out  of  the  original  TE3  dataset  containing  173  P-wave  FM  polarities.  The  analysis  was  repeated 
on  the  smaller  datasets  and  the  results  are  shown  in  Figure  S5  available  in  the  electronic 
supplement  to  Nayak  and  Dreger  (2015).  For  both  data  subsets,  our  inversion  scheme  can 
compute  an  NSS  nearly  equivalent  to  one  estimated  by  the  forwarding  modeling  approach. 
Figure  S5  also  shows  that  ideally,  a  large  quantity  of  FM  polarity  data  is  required  to  reliably 
constrain  source-type  of  any  seismic  event.  However,  it  is  noted  that  the  greatest  benefit  of  our 
approach  is  the  combination  of  waveforms  and  FM  polarities,  which  as  shown  in  previous 
studies  (Ford  et  al.,  2012;  Chiang  et  al.,  2014)  only  a  few  (3  to  10)  FM  polarity  observations  can 
greatly  enhance  the  recovery  of  the  source-type  of  events  like  nuclear  explosions  in  sparse 
monitoring  conditions. 

3.1.7  Joint  Waveform  and  First-Motion  NSS-inversion 

With  the  waveform  and  first-motion  approaches  it  is  now  possible  to  perform  a  joint 
NSS-inversion  of  waveform  and  FM  polarity  data  for  events  TE1,  TE2  and  TE3.  The  two  types 
of  data  are  first  inverse  weighted  by  their  sum-of-squares  to  account  for  the  difference  between 
their  amplitudes.  We  use  100  random  initial  models  for  each  A-spccilic  inversion  and  test  them 
against  the  data.  30  models  with  the  best  fits  are  then  used  as  starting  models  in  the  inversion. 
The  damping  parameter  from  trial  and  error  is  set  to  10.  Figure  7  compares  NSS  computed  using 
wavefonn  and  FM  polarity  data  by  joint  inversion  and  forward  modeling  80  million  MT,  which 
are  very  similar.  For  events  TE1  and  TE2  (Figure  7a, b),  all  FM  polarities  are  positive  and  the 
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inclusion  of  FM  polarities  in  the  NSS  constrains  the  best- fitting  source  type  to  be  dominantly 
volume-increase.  The  maximum  VR  recovered  from  the  NSS  computed  using  the  joint  inversion 
of  wavefonns  and  FM  polarities  for  events  TE1  and  TE2  (84.0%  and  91.9%,  respectively),  are 
very  close  to  the  mean  of  separate  maximum  VR  values  for  FM  polarity  data  (100%  for  both 
events)  and  waveforms  (68.2%  and  85.4%,  respectively)  which  demonstrates  the  success  of  the 
joint  inversion  scheme  proposed  in  this  study.  For  event  TE3  (Figure  7c),  the  maximum  VR 
recovered  from  NSS  computed  using  the  joint  inversion  of  waveforms  and  FM  polarities,  59.4%, 
is  less  than  mean  of  maximum  VR  values  (65.8%)  for  separate  waveform  (80.2%)  and  FM 
polarity  (51.3%)  NSS.  Assuming  NSS  converged  to  the  best-fitting  solution,  this  suggests  that 
MT  solutions  or  source-types  that  best  fit  the  low  frequency  displacement  wavefonns  and  high- 
frequency  P-wave  FM  polarities  separately  are  not  in  perfect  agreement  with  each  other.  While 
FM  polarities  reflect  the  source  mechanism  at  the  beginning  of  the  source-time  history  of  the 
earthquake,  low  frequency  wavefonns  reflect  the  average  source  mechanism  of  the  event  under  a 
point  source  assumption  in  time  and  space,  which  might  lead  to  disagreement  between  best¬ 
fitting  MT  solutions  for  the  two  data-types  in  case  of  complex  events.  For  example,  initial 
radiated  energy  in  explosions  is  expected  to  be  isotropic  and  explosive  leading  to  all  compressive 
FM  polarities  whereas  tectonic  release  during  the  later  part  of  the  time  history  might  lead  to 
deviatoric  components  in  the  point  source  MT  solution  obtained  from  waveform  inversion.  Man¬ 
made  explosions  are  suitable  for  joint  analysis  of  waveforms  and  FM  polarity  data,  because  the 
initial  mechanism  must  be  explosive  and  volume-increase,  leading  to  positive  first  motions 
irrespective  of  the  nature  of  secondary  mechanism  responsible  for  non-isotropic  components  that 
are  commonly  found  in  MT  solutions  of  explosions.  For  events  in  geothennal  and  volcanic 
environments,  source-time  functions  could  be  complex  and  composite  mechanisms  could  be 
diverse,  and  therefore,  it  might  be  difficult  to  reconcile  FM  polarity  data  with  composite  point 
source  MT  solutions.  It  is  therefore  expected  that  joint  inversion  of  polarity  data  and  low- 
frequency  wavefonns  works  best  for  small  events  with  short  duration  impulsive  source-time 
functions.  Nevertheless,  the  analysis  can  identify  potential  discrepancies  between  the  two  data 
sets  thereby  raising  caution  and  a  need  for  further  analysis. 

For  all  events  and  all  types  of  NSS,  the  maximum  VR  as  well  the  size  of  the  region  of 
best-fitting  MT  solutions  (highest  contour  levels  in  Figures  4,  6  and  7)  recovered  by  the  NSS- 
inversion  approach  are  either  same  as  or  greater  than  those  recovered  by  forward  modeling  80 
million  random  MT.  For  the  waveform-only  NSS  the  VR  recovered  by  the  NSS-inversion 
approach  are  the  either  same  as  or  greater  than  those  recovered  by  forward  modeling  random  MT 
for  almost  all  source-types  (Figure  4).  For  the  NSS  computed  using  wavefonn  and  polarity  data 
simultaneously,  we  recognize  there  are  few  eigenvalue  sets  for  which  VR  values  of  the  best¬ 
fitting  MT  solutions  recovered  by  our  inversion  method  are  not  as  good  as  those  from  the  large 
population  of  random  MT  solutions  (for  example,  see  outennost  contours  of  Figure  7b).  These 
differences  are  possibly  caused  by  the  inability  of  the  inversion  to  converge  further  towards  the 
true  solutions  on  account  of  the  derivatives  of  polarity  data  having  near-zero  values  (see  Figure  5 
and  its  explanation  earlier  in  this  section).  Comparing  the  results  of  multiple  runs  of  joint 
waveform  and  FM  polarity  NSS  inversion,  we  have  observed  that  while  the  shape  of  VR 
contours  in  NSS  are  similar,  the  maximum  VR  recovered  can  vary  by  1-2%.  Further  work  is 
required  to  design  an  improved  inversion  scheme  with  a  better  way  of  incorporating  polarities 
that  converges  to  the  true  global  VR  maximum  at  each  instance.  The  computation  of  NSS  by 
joint  inversion  of  wavefonns  and  polarity  data  also  takes  2-4  times  longer  than  the  computation 
when  inverting  waveforms  alone.  However,  this  increase  in  computation  time  can  be  remedied 
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by  not  inverting  at  source-type  grid  points  where  the  waveform-only  VR  is  less  than  a  particular 
threshold  (for  example,  40%). 


£*P'osion  NSS  Inversion 
^  I  Event  TE1 


NSS  Random 
Event  TE1 


+Crack 


+Dipole. 


+CLVD 


■CLVD 


■Crack 


Implosion 


NSS  Random 
Event  TE2 


NSS  Inversion 
Event  TE2 


NSS  Random 
Event  TE3 


NSS  Inversion 
Event  TE3 


Figure  7.  NSS  of  the  three  events  using  both  low-frequency  displacement  waveforms  and  FM 
polarity  data:  (a)  TE1,  (b)  TE2  and  (c)  TE3.  Explanation  of  symbols,  shading  and  contours  is 
same  as  in  Figure  4.  The  color  version  of  this  figure  is  available  only  in  the  electronic  edition. 


Approved  for  public  release;  distribution  is  unlimited. 

19 


3.1.8  Conclusions 


In  this  study,  we  define  elements  of  the  general  MT  in  terms  of:  (1)  its  normalized 
eigenvalues  that  characterize  its  source-type,  (2)  a  seismic  moment  scale  factor  that  scales  the 
nonnalized  eigenvalues  to  appropriate  size  or  scalar  moment,  and  (3)  its  eigenvectors  that 
specify  its  orientation.  We  utilize  this  formulation  to  implement  an  iterative  damped  LS 
inversion  scheme  to  invert  displacement  waveforms  for  best  fitting  eigenvectors  and  the  moment 
scale  factor  for  specific  eigenvalues,  which  results  in  the  best-fitting  MT  solution  for  a  specific 
source-type.  However,  the  expressions  are  general  and  can  be  used  with  other  appropriate 
inversion  techniques  as  well.  Our  technique  is  successfully  demonstrated  by  estimating  best 
fitting  MT  solutions  for  synthetic  data  assuming  various  source  types:  cracks,  explosion  and  DC. 
For  low  frequency  displacement  waveforms  of  the  three  example  events,  we  find  the  NSS 
computed  using  the  inversion  approach  to  be  faster  and  more  accurate  by  VR  ~  0-3%  compared 
to  NSS  computed  by  forward  modeling  80  million  randomly  generated  MT.  To  better  constrain 
the  source-types  of  these  seismic  events,  we  employ  an  approximation  of  the  sign  function  in 
order  to  invert  FM  polarity  data  along  with  displacement  waveforms  using  our  derivative-based 
inversion  scheme.  We  find  that  our  inversion  method  is  more  successful  than  the  random  search 
approach  in  recovering  the  MT  solution  with  the  maximum  VR  as  well  as  the  region  of  best- 
fitting  MT  solutions  or  source-types  for  NSS  computed  using  low-frequency  displacement 
waveforms  and  P-wave  FM  polarities  both  separately  and  jointly.  The  inclusion  of  P- wave  first- 
motion  observations  with  long-period  waveforms  narrows  the  range  of  possible  MT  solutions  in 
source-type  space  leading  to  improved  source-type  discrimination.  It  would  be  straightforward  to 
also  incorporate  body  wave  amplitude  ratios,  and  such  a  combination  of  data  types  could  be 
useful  in  cases  where  individual  data  sets  are  very  sparse,  such  as  for  small  magnitude  induced 
earthquakes. 
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3.2  Analysis  of  the  Effects  of  Vanishing  Traction  for  Shallowly  Buried 
Sources 

A  potential  issue  for  moment  tensor  inversion  of  shallow  seismic  sources  is  that  some 
moment  tensor  components  have  vanishing  amplitudes  at  the  free  surface,  which  can  result  in 
bias  in  the  moment  tensor  solution.  The  effects  of  the  free-surface  on  the  stability  of  the  moment 
tensor  method  becomes  important  as  we  continue  to  investigate  and  improve  the  capabilities  of 
regional  full  moment  tensor  inversion  for  source-type  identification  and  discrimination.  It  is 
important  to  understand  these  free  surface  effects  on  discriminating  shallow  explosive  sources 
for  nuclear  monitoring  purposes.  It  may  also  be  important  in  natural  systems  that  have  shallow 
seismicity  such  as  volcanoes  and  geothermal  systems.  In  this  study,  we  examine  the  effects  of 
the  free  surface  on  the  moment  tensor  via  synthetic  testing,  and  apply  the  moment  tensor  based 
discrimination  method  to  the  HUMMING  ALBATROSS  quarry  blasts.  These  shallow  chemical 
explosions  at  approximately  10  m  depth  and  recorded  up  to  several  kilometers  distance  represent 
rather  severe  source-station  geometry  in  terms  of  vanishing  traction  issues.  We  show  that  the 
method  is  capable  of  recovering  a  predominantly  explosive  source  mechanism,  and  the  combined 
waveform  and  first  motion  method  enables  the  unique  discrimination  of  these  events.  Recovering 
the  correct  yield  using  seismic  moment  estimates  from  moment  tensor  inversion  remains 
challenging  but  we  can  begin  to  put  error  bounds  on  our  moment  estimates  using  the  NSS 
technique. 

A  preliminary  report  of  this  section  was  provided  in  the  final  report  for  FA9453-10-C-0263.  It  is 
included  here  to  report  on  the  final  results  of  the  effort,  which  have  been  submitted  for 
publication  in  the  Bulletin  of  the  Seismological  Society  of  America  (Chiang  et  ah,  2016),  and  at 
the  time  of  writing  this  report  is  accepted  pending  minor  revision. 


3.2.1  Introduction 

Wavefonn  inversion  to  determine  the  seismic  moment  tensor  is  now  a  standard  method 
for  determining  the  source  mechanism  of  natural  and  man-made  seismicity,  and  can  identify,  or 
discriminate  different  types  of  seismic  sources.  Such  source-type  identification  is  important  for 
better  understanding  the  physical  processes  of  explosions  as  well  as  other  man-made  seismicity 
and  earthquakes  from  geothermal  Guilhem  et  ah,  2014),  and  volcanic  environments  (Templeton 
and  Dreger,  2006;  Minson  et  ah,  2007)  and  oil  and  gas  operations  (McNamara  et  ah,  2015).  For 
the  nuclear  explosion  discrimination  problem  the  uncertainty  in  a  given  solution  is  as  important 
as  the  best  fitting  parameters  and  it  is  therefore  necessary  to  fully  understand  and  model  possible 
bias  that  can  result  in  such  inversions.  A  potential  issue  for  shallow  seismic  sources  that  are 
effectively  at  the  free  surface  is  that  as  the  traction  vanishes  the  associated  vertical  dip-slip  (DS) 
Green’s  functions  have  vanishing  amplitudes  (Julian  et  ah,  1998,  Stevens  and  Murphy,  2001), 
which  in  turn  can  result  in  the  indeterminacy  of  the  Mxz  and  Myz  components  of  the  moment 
tensor  and  therefore  bias  in  the  moment  tensor  solution.  The  free-surface  effect  was  noted  in  a 
study  on  fundamental  Love  and  Rayleigh  waves  for  nuclear  explosions  and  associated  tectonic 
release  (Given  and  Mellman,  1986).  The  effects  of  the  free  surface  on  the  stability  of  the  moment 
tensor  becomes  important  as  we  continue  to  investigate  and  improve  the  capabilities  of  regional 
full  waveform  moment  tensor  inversion  for  source-type  identification  and  discrimination.  It  is 
important  to  understand  its  effects  for  discriminating  shallow  explosive  sources  for  nuclear 
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monitoring,  but  could  also  be  important  in  natural  systems  that  have  shallow  seismicity  such  as 
volcanic  environments  and  geothermal  systems,  and  other  manmade  shallow  seismicity  related  to 
anthropogenic  activities  such  as  hydraulic  fracturing  and  mining. 

To  investigate  the  potential  issues  that  could  arise  in  the  estimation  of  moment  tensors  for 
shallow  sources,  we  perfonn  a  series  of  synthetic  tests  to  document  and  understand  the  effects  of 
free  surface  vanishing  traction  (FSVT)  on  the  total  seismic  moment,  isotropic  seismic  moment 
and  the  source  mechanism.  We  evaluate  the  sensitivity  of  the  moment  tensor  solutions  as  a 
function  of  source  depth,  data  quality,  frequency  and  velocity  model.  Based  on  what  we  learn 
from  the  synthetic  studies,  we  applied  the  moment  tensor  method  to  the  HUMMING 
ALBATROSS  quarry  blast  events,  this  is  an  excellent  dataset  in  terms  of  understanding  the 
effects  of  FSVT  and  evaluating  the  discrimination  capabilities  of  the  moment  tensor  method  as  a 
function  of  source  depth  and  frequency.  These  small  chemical  explosions  are  approximately  10 
m  deep  and  are  recorded  at  up  to  several  km  distances.  Therefore  the  data  represent  rather  severe 
source-station  geometry  in  terms  of  vanishing  traction  issues.  We  show  that  the  combined 
method  utilizing  both  complete  seismic  waveforms  and  P-wave  first  motion  polarities  is  able  to 
obtain  robust  full  moment  tensor  solutions  that  are  comprised  predominantly  by  an  isotropic  or 
explosive  component.  However,  unlike  the  synthetic  studies,  yield  estimation  using  the  real 
quarry  blast  moment  tensor  inversion  results  remains  challenging. 


3.2.2  Methods 

The  seismic  moment  tensor  consists  of  nine  force  couples  that  represent  the  equivalent 
body  forces  for  seismic  sources  of  different  geometries  (Jost  and  Herrmann,  1989)  that  due  to 
conservation  of  angular  momentum  reduce  to  six  independent  couples  and  dipoles.  The  data  (e.g. 
displacement  waveforms)  are  represented  by  the  convolution  of  Green’s  functions  for  a  given 
Earth  model,  source  term  and  the  moment  tensor  elements.  We  obtain  the  individual  moment 
tensor  elements  by  inverting  the  3 -component,  complete  wavefonn  data  using  a  time  domain, 
generalized  least  square  inversion  (Minson  and  Dreger,  2008),  and  the  goodness  of  fit  between 
the  data  and  synthetics  is  measured  by  the  variance  reduction  (VR): 


VR  = 


gjWj(di-Si)^\ 
XiWjd,2  ) 


x  100, 


(19) 


where  d  is  the  data,  s  is  the  synthetic  waveforms  and  w  is  the  distance  weighting  at  each  station  i. 
Because  the  data  are  linear  combinations  of  the  Green’s  functions  weighted  by  their  associated 
moment  tensor  elements,  one  major  source  of  error  in  the  inverted  moment  tensor  solution  comes 
from  the  assumed  velocity  model.  A  well-calibrated  velocity  model  is  important  to  obtain  robust 
estimates  of  the  source  parameters.  For  HUMMING  ALBATROSS  we  used  the  ID  velocity 
model  by  Saikia  et  al.  (1990)  computed  from  Rg  wave  dispersion,  and  calculated  the  Green’s 
functions  using  frequency-wavenumber  integration  (Wang  and  Herrmann,  1980;  Herrmann  and 
Wang,  1985;  Hernnann,  2013).  The  inversion  method  also  allows  for  small  time  shifts  between 
the  data  and  Green’s  functions  to  compensate  for  uncertainties  in  origin  time,  location,  and 
velocity  structure.  To  properly  account  for  non-double-couple  radiation  we  used  the  Bowers  and 
Hudson  (1999)  fonnulation  to  calculate  the  moment  magnitude. 

To  assess  the  full  uncertainties  in  the  moment  tensor  inversion  and  avoid  the  need  to 
decide  on  a  particular  moment  tensor  decomposition  scheme,  we  examine  the  moment  tensor 
solution  in  terms  of  the  maximum  fit  surface  called  Network  Sensitivity  Solutions  [NSS]  (Ford  et 
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al.,  2010).  The  NSS  presents  the  goodness  of  fit  between  data  and  a  suite  of  moment  tensor 
solutions  for  a  given  station  configuration,  Earth  model,  and  frequency  band.  We  contoured  the 
populations  of  best-fitting  moment  tensor  solutions  that  have  VR  equal  and  above  a  certain 
threshold.  In  descending  order,  the  contours  map  out  the  solutions  that  have  VRs  of  98%,  95%, 
90%,  80%,  70%,  60%,  and  50%  of  the  maximum  VR  in  the  NSS  population.  From  the  NSS  of  a 
given  event  we  can  determine  whether  or  not  the  best-fitting  full  moment  tensor  solution  from 
the  inversion  is  well  resolved  to  make  useful  interpretations  about  the  source.  The  maximum  fit 
surface  is  parameterized  and  represented  in  the  source-type  space  (Hudson  et  al,  1989;  Tape  and 
Tape,  2012a;  Tape  and  Tape,  2012b;  Vavrycuk,  2015).  In  essence  the  source-type  diagram  is  a 
graphical  representation  of  the  full  moment  tensor  eigenvalues.  The  source-type  diagram  has  two 
key  parameters  y  and  8,  which  are  the  longitude  and  latitude  of  a  lune  on  the  unit  sphere  (Tape 
and  Tape,  2012a;  Tape  and  Tape,  2012b): 


tany  - 
cos(3  - 


-A1+2A2+A3 
V3(Ai -A3)  ’ 

Ai+A2+A3 

J  3(Af+A|+A|) 


5  = 


(20) 

(21) 

(22) 


P  is  the  colatitude  and  A4-3  are  the  eigenvalues  of  the  full  seismic  moment  tensor  where  V  >  %  > 
A.3.  Equation  (20)  measures  the  deviation  from  a  pure  shear  dislocation  and  equation  (21)  and 
(3.2.4)  describes  the  volume  change.  In  this  convention  when  5  =  0  (no  volume  change):  y  =  0 
describes  a  pure  double-couple  (DC)  source  and  ±30  describes  a  pure  compensated-linear- 
vector-dipole  (CLVD)  source;  and  y  =  0  and  5  =  ±90  represents  a  spherical  explosion  (±V)  and 
implosion  (-V),  respectively.  Understanding  the  relative  contributions  of  the  different  moment 
tensor  elements  provides  insights  into  the  complex  source  processes  of  explosions  as  well  as 
other  seismic  events.  This  representation  of  the  seismic  source  has  been  shown  to  result  in 
separate  populations  for  explosions,  underground  cavity  collapse  and  earthquakes  (Ford  et  al., 
2009a;  Ford  et  al.,  2009b;  Dreger  et  al.,  2008;  Ford  et  al.,  2008;  Chiang  et  al.,  2014)  enabling 
discrimination  capability. 

Instead  of  the  original  forward-modeling  approach  by  Ford  et  al.  (2010)  we  implemented 
the  damped  least  square  inversion  scheme  by  Nayak  and  Dreger  (2015)  to  compute  the  NSS. 
This  method  significantly  reduces  the  computation  time,  and  recovers  a  true  maximum  fit  surface 
in  the  source-type  space.  In  addition  to  long  period  wavefonn  data  we  include  P-wave  first 
motion  polarities  in  the  NSS  calculations  as  additional  constraints  (Ford  et  al.,  2012;  Guilhem  at 
al.,  2014;  Chiang  et  al.,  2014).  The  damped  least  square  inversion  finds  the  moment  tensor 
solutions  that  best  fit  both  the  complete  waveform  data  and  the  P-wave  first  motion  polarities. 
Similar  to  comparing  synthetic  waveforms  to  data,  we  compare  the  theoretical  P-wave  first 
motions  against  the  observed  polarities.  We  assign  -1  for  downward  motion  and  +1  for  upward 
motion,  and  each  observed  polarity  is  weighted  by  the  quality  of  the  pick.  The  first  motion  VR  is 
calculated  as: 


VR  =  1 


XWj(Pol0bs  Polsynth)  \ 

5>iPoiobs2  J 


x  100. 


(23) 
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The  final  combined  VR  is  computed  by  weighing  the  two  data  sets  equally.  Incorporating  the 
first  motion  data  proves  to  be  a  powerful  tool  in  reducing  solution  uncertainties.  Constraints  from 
first  motion  can  reduce  uncertainties  due  to  the  theoretical  ISO-CLVD  (Isotropic-CLVD) 
tradeoff  (Ford  et  ah,  2012;  Chiang  et  ah,  2014),  sparse  station  coverage  (Chiang  et  ah,  2014), 
illuminate  possible  complexities  in  earthquake  ruptures  in  geothermal  environments  (Guilhem  et 
ah,  2014),  and  as  illustrated  later  in  this  study,  free  surface  effects  at  shallow  depths.  We  find 
that  the  additional  polarity  constraints  assist  by  uniquely  discriminating  the  events  as 
predominantly  explosive  and  greatly  enhance  the  discriminatory  power  of  the  moment  tensor- 
based  method. 

3.2.3  Free-Surface  Vanishing  Traction 

We  generate  the  ten  fundamental  Green’s  functions  at  a  distance  of  1 00-km  from  the 
source  and  with  source  depths  ranging  from  0.2  to  1.2  km.  We  apply  an  acausal  bandpass 
Butterworth  filter  at  10  to  50  seconds  to  filter  the  Green’s  functions.  As  shown  in  Figure  8a  we 
see  a  strong  source  depth  dependency  on  the  vertical  dip-slip  (DS)  fundamental  Green’s 
functions  associated  with  the  Mxz  and  Myz  elements  for  all  three  components:  vertical  (ZDS), 
radial  (RDS)  and  transverse  (TDS)  in  which  there  is  a  systematic  reduction  in  displacement 
amplitude  with  decreasing  source  depth.  The  averaged  Fourier  spectral  amplitude  shows  a  linear 
relationship  between  the  source  depth  and  DS  wavefonn  amplitudes  (Fig.  8b).  In  contrast,  the 
vertical  strike-slip  Green’s  functions  for  all  three  components  (ZSS,  RSS  and  TSS)  and  the 
explosion  Green’s  functions  for  the  vertical  and  radial  components  (ZEX,  REX)  show  little  to  no 
variation  in  amplitude  and  wavefonn  with  respect  to  decreasing  source  depth.  The  vertical  and 
radial  45-degree  dip-slip  Green’s  functions  (ZDD  and  RDD)  show  variations  in  amplitude  as 
well  due  to  the  constructive  and  destructive  interference  of  waves  interacting  with  the  free 
surface.  While  the  wave  interference  appears  minor  in  the  10  to  50  second  period  passband  (Fig. 
8a)  it  is  more  pronounced  in  the  unfiltered  displacement  Green’s  functions. 
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a.  Theoretical  Green’s  Functions 
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Figure  8.  a)  Green’s  functions  (GF)  computed  using  the  Song  et  al.,  (1996)  ID  western  U.S. 
velocity  model.  GF  are  in  displacement  and  filtered  between  10-50  second  period.  The  traces  for 
each  component,  from  top  to  bottom,  are  at  200,  400,  600,  800,  1000,  and  1200  m  source  depths, 
b)  Averaged  spectral  displacement  amplitudes  between  10  to  50  seconds  for  the  ten  fundamental 
Green’s  functions,  c)  Moment  tensor  elements  computed  using  the  same  velocity  model  and 
filtered  between  10-50  second  period. 
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Since  traction  perpendicular  to  the  vertical  vanishes  at  the  free  surface,  the  inversion  may 
not  resolve  the  Mxz  and  Myz,  as  well  as  the  isotropic  components  Mxx,  Myy,  Mzz  of  the  moment 
tensor.  It  is  important  to  note  however  that  while  there  are  strong  effects  on  amplitude,  the 
waveforms  remain  similar  and  there  is  little  effect  on  the  phase  of  the  waveforms  of  these 
components.  The  systematic  behavior  of  the  Green’s  function  suggests  that  it  may  be  possible  to 
scale  the  isotropic  moment  (Miso)  to  recover  the  correct  moment,  enabling  a  more  robust 
estimate  of  the  explosive  yield.  But  because  the  ability  to  resolve  seismic  waveforms  at  the  free 
surface  depends  on  the  velocity  model,  frequency  band,  station  configuration  and  background 
noise  level,  it  is  necessary  to  understand  the  behavior  of  the  moment  tensor  inversion  for  each 
particular  monitoring  scenario  in  order  to  correct  for  bias  in  the  seismic  moment.  Poor  station 
coverage  and  low  signal-to-noise  ratio  (SNR)  will  decrease  the  method’s  ability  to  resolve  the 
observed  waveforms. 

In  the  synthetic  study  we  examine  the  ability  to  recover  the  correct  moment  tensor 
elements  by  generating  a  suite  of  velocity  models  derived  from  a  ID  western  U.S.  velocity 
reference  model  (Song  et  ah,  1996).  We  choose  this  particular  model  as  a  reference  because  it 
was  used  to  model  the  NTS  explosions  in  Ford  et  al.  (2009).  We  generated  the  suite  of  models  by 
splitting  the  top  2.5-km  thick  layer  in  the  reference  model  into  two  separate  layers.  We 
systematically  adjust  the  thickness  and  velocity  of  the  two  new  layers  (Figure  9),  but  constrain 
the  perturbations  of  the  two  parameters  (velocity  and  layer  thickness)  by  maintaining  the  same 
vertical  travel  time  as  the  reference  model.  The  reason  to  maintain  the  same  vertical  travel  time 
is  to  generate  different  but  comparable  velocity  models  and  to  minimize  travel-time  differences. 
For  each  ID  model,  we  generate  displacement  Green’s  functions  at  regional  distances  between 
100  and  400  km  with  source  depths  ranging  from  0.2  to  3.5  km.  Using  the  same  set  of  Green’s 
functions  we  generate  two  types  of  synthetic  data  with  different  source  mechanisms:  1)  a  pure 
explosion  case  (EXP)  and  2)  a  composite  source  case  (DC  and  EXP)  with  a  DC  to  EXP  ratio  of 
2:3.  We  compute  the  synthetic  data  following  the  expressions  from  Minson  and  Dreger  (2008), 
and  add  random  Gaussian  white  noise.  We  scale  the  amplitude  of  the  random  noise  for  a  given 
SNR  using  the  following  equation: 

RMSSjBnai 

SNR  =  201oglo(- - — S; - ).  (24) 

^noise^'^noise 

The  scaling  factor  Cnoise  is  calculated  for  a  given  SNR  using  the  root  mean  squared  (RMS) 
amplitudes  of  the  signal  and  noise,  then  the  noise-included  data  and  Green’s  functions  are 
bandpass-filtered  between  10  to  50  seconds  using  an  acausal  Butterworth  filter.  To  model  after 
real-life  monitoring  capabilities  we  use  an  SNR  of  10  but  also  tested  different  SNRs.  We 
implemented  a  semi-ideal  four-station  coverage  for  the  moment  tensor  inversion,  consisting  of 
source-to-receiver  distances  at  100,  200,  300  and  400  km,  and  distributed  in  semi-regular 
azimuths. 
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Figure  9.  Velocity  models  derived  from  the  Song  et  al.  (1996)  ID  model  by  keeping  the  top  2.5- 
km  vertical  travel  time  constant.  The  model  parameters  are  the  same  below  2.5-km  depth. 

In  this  test  we  keep  the  source  depths  of  the  synthetic  data  and  the  Green’s  functions  used 
in  the  inversion  the  same  to  isolate  the  effect  of  only  vanishing  traction  for  shallow  depth  of 
burial.  Of  the  59  velocity  models  tested  the  full  moment  tensor  inversion  successfully  recovers 
the  correct  mechanism  for  both  the  pure  explosion  case  and  the  composite  case  over  the  targeted 
depth  range  (<  1  km)  for  nuclear  explosions,  as  well  as  at  deeper  depths.  Given  the  same  filter 
parameters  and  station  configuration,  models  that  have  a  thick  continuous  layer  generally  have 
larger  deviations  in  the  moment  estimates,  whereas  the  more  gradient-like  models  show  little  to 
no  change  in  their  moment  estimates  (Figure  10).  At  the  shallowest  depths  (<  0.5  km),  for  the 
pure  explosion  case  the  Miso  estimates  all  fall  within  about  10%  of  the  input  values  (Figure  10a), 
whereas  the  total  moment  varies  even  less.  The  composite  case  exhibits  greater  deviations  from 
the  input  value  at  the  shallowest  source  depths,  however  all  moment  estimates  are  within  about 
20%  of  the  input  values  (Figure  10b).  In  general  the  seismic  moment  estimates  approach  the 
correct  value  as  the  source  depth  increases.  The  IV1ZZ  component  is  not  as  well  constrained,  thus 
affecting  the  method’s  ability  to  resolve  the  Miso-  In  most  cases  the  FSVT  has  little  effect  on 
recovering  the  correct  mechanism  for  models  with  a  shallow  velocity  gradient.  However,  as  the 
noise  level  increases  the  bias  in  the  seismic  moment  and  mechanism  also  increases.  The  key 
observation  from  the  synthetic  tests  is  that  the  bias  in  the  recovered  moment  and  source 
mechanism  of  very  shallow  sources  depends  on  various  factors,  including  the  velocity  model. 
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Figure  10.  Isotropic  moment  and  total  seismic  moment  percent  change  for  a)  a  pure  explosion 
source  and  b)  a  composite  source,  plotted  as  a  function  of  source  depth.  The  average  value  of  all 
59  models  is  the  solid  red  line,  the  shaded  region  is  2-sigma  from  the  mean,  and  the  dashed  line 
represents  no  deviation  from  the  input  seismic  moments. 


In  practice  when  we  perform  moment  tensor  inversions  the  source  depth  is  also  tested 
because  we  do  not  know  the  true  source  depth.  Thus  we  also  looked  at  synthetic  comparisons 
where  we  fixed  the  source  depths  to  0.2,  0.5,  1.5  and  3.0  km  and  for  each  depth  the  full  suite  of 
Green’s  function  depths  are  tested.  The  results  for  both  source  types  are  shown  in  Figure  1 1.  In 
this  synthetic  test  the  deviation  from  the  true  seismic  moment  is  much  greater  due  to  errors  from 
incorrect  source  depth,  the  percent  change  in  total  and  isotropic  moment  decrease  as  the  Green’s 
function  depths  approach  the  correct  value.  For  the  four  source  depths  we  examined,  the  mean 
percent  changes  in  moments  are  similar  but  the  spreads  are  different.  For  example  when  source 
depth  is  at  1 .5  km  the  change  in  isotropic  moment  away  from  the  true  source  depth  encompasses 
a  much  wider  range  compare  to  other  tested  source  depths.  The  change  in  VR  is  most  prominent 
at  source  depth  of  3.0  km  because  a  boundary  layer  exists  at  2.5  km  for  all  the  velocity  models. 


Approved  for  public  release;  distribution  is  unlimited. 
28 


a.  Pure  Explosion 
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Figure  1 1 .  Fixed  Source  Depth  Sensitivity  Analysis.  Source  depth  is  fixed  while  we  tested  the 
full  suite  of  Green’s  function  depths  (x-axis),  and  compare  the  VR,  total  moment  percent  change 
and  isotropic  moment  percent  change  as  a  function  of  varying  Green’s  function  depth  for  (a)  a 
pure  explosion  source  and  (b)  a  composite  source.  Similarly,  the  solid  lines  are  the  average 
values  of  all  59  models  and  the  shaded  regions  are  2-sigma  from  the  mean.  The  lines  and  shaded 
regions  are  color-coded  according  to  source  depths  (data). 


3.2.4  HUMMING  ALBATROSS 

The  synthetic  tests  show  although  the  loss  of  traction  at  the  free  surface  can  affect  the 
method’s  ability  to  resolve  the  moment  at  shallow  depths,  the  impact  on  the  recoverability  of  the 
source  mechanism  is  much  less.  We  evaluate  the  performance  of  the  moment  tensor  method 
using  the  data  from  the  HUMMING  ALBATROSS  series.  The  data  sets  consist  of  twenty-six 
broadband,  strong  motion  and  short-period  seismic  recordings  (Figure  12).  The  chemical 
explosions  are  denoted  at  very  shallow  depths  (<  20  m)  and  recorded  at  stations  up  to  the  several 
kilometers.  We  apply  the  moment  tensor  based  discrimination  method  to  three  of  the  five 
chemical  explosions,  the  two  smallest  events  in  the  series  have  poor  SNR  below  3  Hz  and  we 
cannot  estimate  the  source  parameters  with  high  confidence.  Weston  Geophysical  Corp. 
provided  the  instrument  corrected  velocity  and  acceleration  data  for  HUMMING  ALBATROSS, 
and  we  prepare  the  velocity  waveforms  for  source  analysis  by  rotating  the  horizontal  recordings 
to  radial  and  transverse  components,  and  bandpass-filtered  the  data  and  Green’s  functions  with 
an  acausal  Butterworth  filter.  Depending  on  the  instrument  response  and  data  quality  the 
bandpass  filter  applied  to  each  station  ranges  between  0.4  to  1.25  seconds.  Table  2  lists  the 
source  depth  and  filter  corners  used,  and  the  moment  magnitude  (Mw),  MiSo  and  wavefonn  fits 
(VR)  from  time  domain  full  moment  tensor  inversions. 
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Figure  12.  Event  and  station  locations  for  the  HUMMING  ALBATROSS  series.  The  seismic 
array  includes  broadband  instruments  (square),  short  period  sensors  (triangle)  and  accelerometers 
(inverted  triangle).  We  looked  at  three  of  the  five  chemical  explosions  (circle)  and  the 
background  colors  represent  the  local  topography  where  green  is  lower  elevation. 


Table  2.  Moment  Tensor  Solution  and  Yield 


Event 

Centroid 

Depth 

(m) 

Filter  (s) 

Mw 

MIS0  (dyne- 
cm) 

VR 

(%) 

Estimated 

Yield 

(tons) 

Yield 

(tons) 

Production 

Varied 

0.5-1.25 

2.5 

3.24747  x  10iy 

70 

28.943 

10.189 

800 

9 

0.5-0.83 

1.9 

3.36495  x  1018 

72 

2.999 

0.394 

400 

11 

o 

1 

1.4 

9.91028  x  101V 

64 

0.883 

0.201 

Following  the  results  of  Ford  et  al.  (2012)  and  Chiang  et  al.  (2014)  we  utilize  full 
wavefonn  data  from  broadband  and/or  short  period  stations  for  the  time  domain  wavefonn 
inversion,  and  include  both  waveform  and  P-wave  first  motion  polarities  from  all  twenty-six 
stations  into  our  source-type  uncertainty  analysis.  We  present  results  from  moment  tensor 
sensitivity  analysis  both  as  a  function  of  depth  and  frequency  for  the  largest,  delayed-fire 
production  shot. 

The  time-domain  moment  tensor  inversion  yield  similar  results  for  all  three  explosions  in 
which  the  full  moment  tensor  solutions  all  have  relatively  large  isotropic  components  that  are 
approximately  between  40  to  50  percent  of  the  total  seismic  moment,  and  the  deviatoric 
inversions  all  result  in  vertical  dip-slip-like  mechanisms  (Figure  13).  The  deviatoric  solutions  are 
thus  dominated  by  the  Mxz  and  Myz  components,  which  as  shown  in  Figure  8c,  the  TDS,  RDS 
and  ZDS  Green’s  functions,  show  a  marked  reduction  in  amplitude  with  decreasing  source  depth. 
An  example  of  the  waveform  fits  and  moment  tensor  solutions  for  the  production  shot  is  shown 
in  Figure  14.  All  three  explosions  have  strong  shear  wave  energy  in  the  transverse  component, 
and  many  of  them  have  amplitudes  comparable  to  Rayleigh  waves  in  the  vertical  component.  We 
avoid  using  stations  in  the  direction  of  the  free  face  (quarry  cliff)  for  the  production  shot  because 
of  spallation  at  the  cliff  face.  For  shot  800  we  tried  a  combination  of  stations  with  and  without 
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those  in  the  direction  of  the  free  face  and  found  the  moment  tensor  solutions  to  be  very  similar. 


Figure  13.  HUMMING  ALBATROSS  moment  tensor  solutions  from  time  domain  full  and 
deviatoric  inversion,  and  best  solution  from  combined  wavefonn  and  P-wave  first  motion  (filled 
circle)  Network  Sensitivity  Solutions.  All  polarities  are  up. 


The  production  shot  source  parameter  uncertainties  using  only  seismic  wavefonns  show 
the  best-fitting  mechanisms  are  mostly  isotropic  (Figure  15a).  The  source  parameters  are  well 
constrained  for  solutions  with  VRs  equal  and  above  95%  of  the  maximum  VR,  but  if  we  relax 
the  threshold  down  to  90%  of  the  maximum  VR  the  population  of  well-fitted  solutions  extends 
across  the  DC/deviatoric  line.  Previous  NSS  studies  of  nuclear  explosions  have  shown  tradeoffs 
between  ISO  (+V)  and  -CLVD  mechanisms,  but  generally  not  between  ISO  and  DC  mechanisms 
(Ford  et  ah,  2012;  Chiang  et  ah,  2014).  A  combination  of  strong  Love  waves  and  free-surface 
effects  most  likely  contributed  to  this  apparent  ISO-DC  tradeoff  and  it  is  especially  pronounced 
for  the  two  smaller  shots.  To  increase  the  confidence  in  our  moment  tensor  solutions  we  bring  in 
constraints  from  P-wave  first  motion  polarities.  We  have  excellent  azimuthal  coverage  for  the 
first  motion  observations  that  sampled  the  entire  focal  sphere,  significantly  reducing  the  NSS  to  a 
predominantly  explosive  mechanism  (Figure  15b).  The  best  solutions  from  the  combined  NSS 
are  similar  to  the  solutions  from  the  waveform-only  inversions  (Figure  13)  Moment  tensor 
solutions  that  deviate  away  from  a  theoretical  opening  crack  do  not  fit  the  observed  polarities. 
NSSs  for  shots  800  and  400  show  similar  behaviors,  namely  the  ISO-DC  tradeoff  using  only 
waveforms,  and  a  well-constrained  solution  when  combining  both  waveform  (dominated  by 
surface  waves)  and  P-wave  first  motion  polarities.  We  also  present  an  example  of  a  first  motion 
only  NSS  (Figure  15c),  the  maximum  fit  surface  using  only  polarity  data  is  extensive,  therefore 
first  motions  only  are  not  sufficient  for  event  source-type  discrimination. 
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Figure  14.  Production  shot  focal  mechanisms  and  waveform  fits  from  full  and  deviatoric  moment 
tensor  inversions. 


Given  the  frequency  band  and  source-receiver  distance,  the  moment  tensor  inversion  is 
not  very  sensitive  to  source  depths  shallower  than  ~100  m.  We  estimate  source  depth  sensitivity 
by  looking  at  the  VR  as  a  function  of  source  depth  and  observe  that  the  deviatoric  VR  drops 
sharply  starting  around  100  m,  whereas  the  full  VR  shows  a  more  gradual  decrease  in  VR.  If  we 
have  no  prior  knowledge  on  source  depth,  the  sensitivity  analysis  shows  we  can  constrain  the 
depth  to  be  shallower  than  - 1 00  m,  indicative  of  possible  manmade  seismicity  since  natural 
earthquakes  rarely  occur  at  these  depths.  Although  we  do  not  have  much  sensitivity  at  the  very 
shallow  depths,  the  mechanisms  remain  stable  and  predominantly  explosive  at  the  borehole 
depths  (Figure  16a-b).  We  start  to  see  a  tradeoff  between  the  mechanism  and  incorrect  source 
depth  at  around  100  m  for  the  production  shot,  and  around  50  m  for  800  and  400;  however  the 
combined  waveform  and  first  motion  inversion  eliminates  this  tradeoff.  Because  we  know  the 
depth  of  the  borehole  where  the  explosions  are  detonated,  we  focus  on  the  behavior  of  seismic 
moment  at  depths  less  than  and  equal  to  20  m.  Due  to  free  surface  effects  the  total  and  isotropic 
moment  increases  as  source  depth  decreases  (Figure  16c-d),  and  the  increases  in  seismic 
moments  are  controlled  by  changes  in  the  Mzz  and  Myz  components  (Figure  16d),  components 
that  we  anticipated  to  be  strongly  influenced  by  the  free  surface.  The  behavior  of  the  moment 
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tensor  solution  as  a  function  of  depth  is  similar  for  shots  800  and  400,  although  the  free  surface 
effects  are  more  pronounced  at  the  shallowest  depths,  where  the  full  moment  tensor  solutions 
become  vertical  dip-slip  due  to  a  greater  increase  in  moment  from  the  Myz  component. 


Figure  15.  Network  Sensitivity  Solutions  (NSS)  using  a)  waveforms,  b)  waveform  and  first 
motions,  and  c)  first  motions.  Best  MT  solutions  from  time-domain  inversion  (plus)  and  NSS 
(circle)  are  plotted  in  a)  and  b).  Production  shot  NSS  are  plotted  as  shaded  contours,  and 
highlighted  contours  are  populations  of  solutions  with  normalized  VR  of  95%  relative  to  the 
maximum. 

The  Earth  becomes  very  heterogeneous  at  high  frequencies  where  the  smoothed  ID 
model  is  no  longer  adequate  in  characterizing  the  complexities  in  high  frequency  wave 
propagation  between  the  source  and  receiver.  Therefore  the  reason  for  utilizing  low  frequency 
wavefonn  data  is  to  minimize  the  bias  from  unaccounted  velocity  structures.  However,  the 
challenge  for  moment  tensor  inversion  at  these  very  shallow  depths  for  the  HUMMING 
ALBATROSS  data  is  that  the  free  surface  effect  not  only  increases  with  decreasing  source  depth 
but  also  at  longer  periods.  The  production  shot  has  relatively  good  SNR  at  periods  up  to  2.5 
seconds,  but  if  we  fix  the  source  depth  to  11  m  and  look  at  the  free  surface  effects  on  the 
inversion  as  a  function  of  filter  passband  we  see  the  solutions  become  vertical  dip-slip  at  longer 
periods  (Figure  17a-b).  The  total  moment  increases  (Figure  17c)  as  we  go  towards  longer  periods 
due  to  rapidly  increasing  contributions  from  the  Mxz  and  Myz  components  (Figure  17d).  Hence 
finding  the  optimal  filter  passband  that  minimizes  both  the  effects  of  the  velocity  structure  at 
high  frequencies,  and  the  free  surface  at  low  frequencies  is  especially  important  and  potentially 
challenging  for  moment  tensor  analysis  of  shallow  sources. 
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Figure  16.  Production  shot  moment  tensor  solutions  as  a  function  source  depth,  using  only 
waveform  data,  a)  VR  from  full  moment  tensor  inversion  (focal  mechanism)  and  deviatoric 
moment  tensor  inversion  (cross),  b)  full  moment  tensor  decomposition.  Shaded  bars  are  the 
percent  DC,  CLVD  and  ISO  components,  c)  total  moment  from  full  inversion  (square)  and 
deviatoric  inversion  (cross),  and  d)  moment  tensor  elements  (filled  circles)  and  isotropic  moment 
(open  circle). 
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Figure  17.  Production  shot  full  moment  tensor  solutions  as  a  function  of  filter  passband,  using 
only  waveform  data,  a)  VR  and  focal  mechanism),  b)  decomposition,  c)  total  moment  (square) 
and  isotropic  moment  (cross),  and  d)  moment  tensor  elements  (filled  circle).  Horizontal  lines  in 
a)-c)  specify  the  filter  passband  used  for  each  inversion  at  1 1  m  depth. 
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3.2.5  Yield  Estimates 


Empirically  the  magnitude-yield  relationship  can  be  described  by  the  following  equation: 


M  =  A  log10  Y  +  C, 


(25) 


where  M  is  the  seismic  magnitude,  Y  is  the  yield  in  kilotons  and  A  and  C  are  constants  that 
depend  on  the  magnitude  measurement  used.  C  may  also  include  additional  corrections  to 
account  for  differences  in  source  medium  that  are  independent  of  the  yield.  Empirical 
relationships  to  estimate  yield  have  been  developed  using  body  wave  magnitudes  (mb)  and 
surface  wave  magnitudes  [Ms]  (e.g.  Mueller  and  Murphy  1971;  Murphy,  1977;  Burger  et  ah, 
1987;  Murphy,  1996)  from  past  nuclear  explosions  with  known  yields.  Given  and  Mellman, 
(1986)  and  Ekstrom  and  Richards,  (1994)  performed  moment  tensor  inversions  on  the  Shagan 
River  explosions  from  the  fonner  Soviet  Union  and  explosions  from  the  U.S.  Nevada  Test  Site 
(NTS),  and  computed  the  moment-yield  relationships.  Yield  estimation  using  the  Miso  from 
moment  tensor  inversions  can  be  challenging  due  to  difficulties  assessing  the  contributions  from 
the  isotropic  and  the  non-isotropic  radiation  of  the  source  as  discussed  in  Stevens  and  Murphy 
(2001),  as  well  as  due  to  free  surface  effects  at  shallow  depths.  However  as  demonstrated  in  Ford 
et  al.  (2012)  and  Chiang  et  al.  (2014)  and  the  previous  sections,  source-type  analysis  (NSS),  and 
the  combination  of  waveform  and  first  motion  data  can  help  to  better  constrain  the  isotropic 
radiation  of  the  moment  tensor. 

We  calculated  the  yields  for  HUMMING  ALBATROSS  using  the  waveform  moment 
tensor  results  at  the  reported  centroid  depths,  and  the  Miso-yield  relationship  Miso  =  logioY  + 
14.05  from  Stevens  and  Murphy  (2001).  The  empirical  relationship  is  derived  from  previously 
published  moment  tensor  inversion  results,  and  announced  yields  or  estimated  yields  from 
Shagan  and  NTS  explosions  (Fig.  8).  For  comparison,  we  also  calculated  the  linear  regression  by 
solving  both  the  slope  and  intercept  (A  and  C  in  eq.  3.2.7),  using  the  same  explosion  dataset.  To 
use  Stevens  and  Murphy  (2001)’s  empirical  relationship  we  extrapolate  the  curve  down  to  very 
low  values  of  yield  and  isotropic  moment.  Because  the  empirical  relationship  is  derived  from 
nuclear  explosions,  differences  in  the  seismic  coupling  efficacy  between  nuclear  and  chemical 
explosions  (Murphy,  1996)  can  also  bias  our  yield  estimates.  We  present  our  yield  estimates  and 
the  true  yields  in  Table  2.  The  moment  tensor  inversion  at  shallow  depth  overestimates  the  yields 
for  all  three  explosions;  the  estimated  values  are  2.8  times,  7.6  times  and  4.4  times  greater  than 
the  true  yield  for  the  production  shot,  800,  and  400,  respectively.  However,  from  our  combined 
NSS  analysis,  we  can  estimate  the  errors  associated  with  our  moment  tensor  solutions.  We 
choose  a  population  of  moment  tensor  solutions  with  normalized  VRs  >  90%  and  computed  the 
mean  and  2-sigma  from  the  mean.  For  the  Production  shot,  the  mean  is  very  close  to  the  Stevens 
and  Murphy  (2001)  scaling,  and  for  the  two  smaller  shots  the  lower  error  bounds  overlaps  with 
the  empirical  scaling. 

From  our  synthetic  studies  on  the  effects  of  free  surface  at  regional  distances,  we  observe 
the  recovered  mechanism  is  less  sensitive  to  the  source  depth,  but  there  can  be  bias  in  the 
moment  due  to  the  loss  of  traction  at  the  surface.  Although  the  deviations  from  the  true  moment 
are  all  within  20%  from  our  synthetic  tests,  the  amount  of  deviation  also  depends  on  station 
coverage,  frequency  passband  and  station  quality.  Additional  synthetic  tests  on  varying  the 


Approved  for  public  release;  distribution  is  unlimited. 

35 


station  coverage,  frequency  band  and  noise  level  suggest  that  although  the  amount  of  bias  in 
Miso  varies,  the  asymptotic  behavior  as  a  function  of  depth  remains  similar,  and  that  the 
recovered  Miso  approaches  the  true  value  as  the  source  depth  increases.  Thus  in  order  to 
minimize  the  vanishing  traction  effects,  it  may  be  possible  to  invert  for  the  source  parameters  at 
a  deeper  depth  where  the  synthetic  tests  show  the  Miso  approaches  the  true  value.  Here  we 
investigate  the  possibility  to  recover  the  correct  yield  by  intentionally  inverting  for  the  source  at 
a  deeper  depth. 

We  first  perform  similar  synthetic  tests  for  a  pure  explosion  scenario  using  the  production 
shot  station  configuration,  filter  passband  and  ID  velocity  model  to  predict  the  depth  range  in 
which  we  can  minimize  both  the  effects  from  vanishing  traction  and  incorrect  source  depth.  We 
fix  the  depth  of  the  data  but  vary  the  depth  of  the  Green’s  functions  (as  what  is  done  in  practice 
when  source  depth  is  unknown  or  has  large  uncertainties),  and  found  that  the  moment  tensor 
inversion  is  able  to  recover  the  correct  focal  mechanism  and  seismic  moment  at  depths  between 
0.02  and  0.05  km.  This  is  an  idealized  case  with  no  time  shift  in  the  inversion,  high  SNR  and 
prefect  knowledge  of  the  Earth  structure,  the  question  remains  whether  this  adjustment  in  the 
source  depth  can  correctly  estimate  the  yield  from  real  explosions  where  we  have  imperfect 
knowledge  of  the  Earth  structure. 

In  Figure  18b  the  yield  estimates  for  HUMMING  ALBATORSS  as  a  function  of  source 
depth  are  compared.  In  general  because  Miso  decreases  as  source  depth  increases,  the  estimated 
yield  also  decreases.  Yields  computed  from  Miso  obtained  at  0.02  and  0.05  km  depth  do  agree 
with  the  true  yield,  but  as  noted  previously  we  start  to  see  tradeoff  between  incorrect  source 
depth  and  source  mechanism  around  50  to  100  km  where  the  moment  tensors  are  no  longer 
predominately  explosive  using  only  wavefonn  data.  Unlike  the  synthetic  case,  with  real  data, 
introducing  additional  errors  from  incorrect  source  depth  can  have  a  greater  impact  on  the 
recoverability  of  the  mechanism,  particularly  for  the  smaller  shots  800  and  400  where  we  use 
higher  frequency  wavefonns. 
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a.  Moment- Yield  Relationship  from  Stevens  and  Murphy  (2001) 
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b.  HUMMING  ALBATROSS 
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Figure  18.  a)  Isotropic  moment  vs.  yield  from  nuclear  explosions  listed  in  Stevens  and  Murphy 
(2001)  and  HUMMING  ALBATROSS  chemical  explosions.  Two  moment-yield  scaling 
relations  are  shown,  one  with  a  slope  of  one  (Stevens  and  Murphy,  2001)  and  one  where  we 
solve  for  both  the  slope  and  intercept.  The  shaded  regions  are  the  95%  confidence  bounds.  The 
mean  isotropic  moment  and  error  bars  (2-sigma)  for  the  three  HUMMING  ALBATROSS  events 
are  computed  from  a  population  of  moment  tensors  with  normalized  VRs  >  90%  in  the  combined 
NSS,  as  shown  in  Figure  15b.  b)  Isotropic  moment  vs.  yield  for  HUMMING  ALBATROSS  at 
various  depths.  The  yields  are  computed  using  the  empirical  scaling  relation  by  Stevens  and 
Murphy  (2001),  and  the  isotropic  moment  (at  the  reported  centroid  depth)  from  waveform 
moment  tensor  inversion  vs.  the  true  yield  is  also  plotted. 


3.2.6  Discussion 

Synthetic  moment  tensor  inversion  tests  show  that  although  the  loss  of  traction  at  the  free 
surface  affects  the  moment  tensor  inversion’s  ability  to  resolve  the  seismic  moment  at  shallow 
depths,  the  impact  on  the  recoverability  of  the  source-type  (explosion,  double-couple,  composite 
mechanism)  is  substantially  less.  This  is  illustrated  by  the  successful  discrimination  of  the  NTS, 
DPRK,  Soviet  JVE,  Chinese  Lop  Nor  and  HUMMING  ALBATROSS  events  as  explosions. 
However,  the  degree  in  which  free-surface  vanishing  traction  affects  the  seismic  moment 
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depends  not  only  on  the  station  geometry,  noise  levels  and  filter  passband  but  also  on  the 
velocity  model.  Free  surface  effects  are  more  pronounced  when  modeling  longer  wavelengths. 

Non-isotropic  radiation  in  the  HUMMING  ALBATROSS  moment  tensor  solution  is  a 
combination  of  both  the  free-surface  effects  and  strong  SH  waves.  Love  wave  generation  from 
chemical  explosions  is  not  an  uncommon  observation,  and  a  number  of  studies  on  quarry  blasts 
suggest  that  the  non-isotropic  radiation  results  from  spall  of  material  from  the  quarry  face 
(Goforth  and  Bonner,  1995;  Bonner  et  al.,  1996),  and  the  interaction  of  the  wavefield  with  the 
severe  topographic  features  of  quarry  benches  (Barker  et  al.,  1993).  McLaughlin  et  al.,  (2004) 
modeled  observed  short  period  Love  and  Rayleigh  waves  from  a  Texas  quarry  and  found  that  the 
horizontal  throw  of  material  from  the  free  face  (spall  model)  is  the  dominant  mechanism.  Our 
moment  tensor  analysis  does  not  take  into  account  the  horizontal  single-force  component  due  to 
spall,  the  collapse  of  the  quarry  cliff  face,  and  for  the  production  shot,  the  temporal  and  spatial 
spread  from  the  delay-firing  scheme.  These  complications  contribute  to  difficulty  in  unique 
source-type  identification  using  long-period  waveforms  only  (Figure  15a).  The  non-isotropic 
component  of  the  inversion  results  is  manifested  as  a  reverse  mechanism  with  either  a  near¬ 
vertical  or  near-horizontal  dipping  layer,  and  a  vertical  dip-slip  mechanism  in  the  deviatoric 
solution.  Because  the  high  frequency  P  waves  may  be  more  sensitive  to  radiated  energy  from  the 
initial  rupture  (Pearce  and  Rogers,  1989;  Guilhem  et  al.,  2014),  including  P-wave  first  motion 
polarities  are  particularly  useful  to  constrain  the  isotropic  component  of  explosive  events  where 
the  first  arrivals  come  from  the  detonated  explosion.  As  illustrated  in  Figure  15,  the  added 
constraints  from  first  motion  observations  greatly  reduce  the  ambiguity  in  source  mechanism 
from  wavefonns  and  result  in  a  unique  identification  of  the  events  as  having  substantial 
volumetric  components  consistent  with  explosions. 

Moment  tensor  based  event  discrimination  for  shallow  sources  is  reliable  when  SNR  is 
high,  especially  with  the  combination  of  data  from  waveforms  and  first  motions.  When  we 
include  errors  in  the  isotropic  moment  estimated  from  the  NSS  analysis,  we  see  although  there  is 
still  a  bias  towards  lager  moment  (800  and  400)  even  at  such  low  yields  and  low  magnitudes,  the 
isotropic  moments  intersect  with  the  empirical  scaling  relation.  Although  source  inversions  at 
greater  depth  reduce  the  effects  of  vanishing  traction,  in  cases  of  extreme  free-surface  effects  and 
complexities  in  the  source  such  as  the  chemical  explosions  from  this  study,  added  errors  from 
incorrect  source  depth  can  impact  event  discrimination  (in  the  case  of  waveform  only  inversion). 
In  addition  to  the  uncertainties  associated  with  the  moment  estimates,  another  factor  that  needs  to 
be  taken  into  consideration  when  estimating  yield  using  magnitude-yield  scaling  laws  is  the 
nature  of  the  explosive  materials.  For  the  same  yield,  nuclear  explosions,  single-fired  buried 
chemical  explosions,  and  ripple-fired  quarry  blasts  exhibit  differences  in  radiated  seismic  energy 
(Murphy,  1996),  therefore  additional  scaling  is  required  to  factor  in  the  differences  in  seismic 
coupling  for  different  types  of  explosions.  All  three  shots  are  partially  confined  but  the 
production  shot  uses  a  ripple-firing  scheme  and  the  two  smaller  shots  are  single-fired  chemical 
explosions.  Traditionally  single-fired  shots  are  designed  to  maximize  the  strength  of  the  seismic 
signal  (Khalturin  et  al.,  1998). 

3.2.7  Conclusions 

Theoretical  vertical  dip-slip  Green’s  functions  associated  with  the  Mxz  and  Myz 
components  have  vanishing  amplitudes  at  shallow  depths  (<  1  km)  due  to  the  loss  of  traction  at 
the  free  surface,  while  the  Green’s  function  waveforms  look  similar  with  little  phase  distortion. 
However  synthetic  calculations  at  shallow  depths  show  moment  tensor  based  event 
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discrimination  is  reliable  and  the  resolvability  of  the  moment  tensor  solution  depends  on  station 
configuration,  noise  level,  the  frequency  band  and  the  velocity  model. 

We  are  able  to  recover  a  predominantly  explosive  source  mechanism  for  the  three 
HUMMING  ALBATROSS  chemical  explosions  from  moment  tensor  inversions.  Although  the 
source-type  uncertainty  analysis  shows  that  we  cannot  uniquely  characterize  the  events  as 
predominantly  explosive  using  only  wavefonn  data,  the  combined  wavefonn  and  first  motion 
method  enables  the  unique  discrimination  of  these  events.  This  method  has  been  applied  in 
previous  studies  which  have  shown  the  inclusion  of  P-wave  first  motions  in  addition  to  full 
waveform  data  eliminates  the  common  ISO-CLVD  tradeoff  (Ford  et  ah,  2012;  Chiang  et  ah, 
2014)  and  reduces  the  uncertainties  of  sparsely  recorded  underground  explosions  with  strong 
Love  waves  and  reversed  Rayleigh  waves  (Chiang  et  ah,  2014).  In  this  study  we  further 
demonstrate  that  incorporating  the  two  data  sets  is  particularly  useful  in  constraining  the 
isotropic  component  of  explosions,  and  the  method  not  only  applies  to  large  events,  but  also 
small  magnitude,  very  shallow  explosions  that  are  effectively  at  the  free  surface.  The 
combination  of  both  low  frequency  full  wavefonn  data  and  high  frequency  P-wave  polarities 
greatly  enhances  the  capabilities  of  the  moment  tensor  source-type  discrimination  method  in 
cases  of  sparse  station  coverage,  strong  Love  waves  and  free-surface  effects. 

The  moment  tensor  method  is  capable  of  event  discrimination,  although  yield  estimation 
using  the  recovered  absolute  seismic  moment  from  moment  tensor  inversion  remains  challenging 
and  can  have  large  uncertainties,  we  can  begin  to  put  error  bounds  on  our  moment  estimates,  and 
therefore  yield,  using  the  NSS  technique  and  combining  waveform  and  first  motion.  Pure 
explosion  synthetic  tests  suggest  source  inversions  at  deeper  depths  reduce  the  free  surface 
effects  on  the  moment,  but  we  cannot  draw  definite  conclusions  using  the  results  from  the 
HUMMING  ALBATROSS  data  set  due  to  not  only  free-surface  effects,  but  uncertainties 
associated  with  imperfect  knowledge  of  the  Earth  structure,  unaccounted  non-isotropic  radiation 
due  to  the  mass  movement  of  the  quarry  face  and  differences  in  seismic  coupling  between 
different  types  of  explosions.  In  addition,  the  moment-yield  relationship  derived  from  larger 
nuclear  explosions  may  not  be  adequate  for  such  low  magnitude  chemical  explosions. 
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3.3.  Synthetic  Analysis  of  3D  Velocity  Structure  Effects  on  Moment  Tensor 
Inversion 

3.3.1  Introduction 

We  have  established  in  the  previous  two  chapters  that  the  regional  seismic  moment 
tensor  (MT)  method  is  a  robust  and  reliable  technique  to  identify  different  seismic  source 
types  under  sparse  monitoring  situations,  at  shallow  source  depths,  and  when  the 
complexities  in  isotropic  source  processes  are  often  not  well  understood  (e.g.  Masse,  1981). 
For  this  study  we  will  explore  possible  biases  associated  with  one  of  the  most  fundamental 
assumptions  in  seismic  source  analysis,  which  is  the  assumption  that  a  ID  velocity  model  is 
a  good  approximation  to  real  Earth  structure. 

The  ground  motion  (u)  recorded  at  the  receiver  is  a  linear  combination  of  the 
source  (seismic  moment  tensor,  M)  and  the  Green’s  tensor  (G),  where  u  =  GM.  The 
Green’s  tensor  represents  the  impulse  response  of  the  medium,  describing  the  seismic 
wave  propagation  between  the  source  and  receiver  for  a  reference  Earth  model.  Basic 
source  types  are  typically  comprised  of  fundamental  fault  representations  (e.g.  Jost  and 
Herrmann,  1989),  but  they  can  also  be  fundamental  single-couples  and  vector  dipoles.  In 
previous  studies  (e.g.  Ford  et  al.,  2009b;  Walter  et  al.,  2009;  Nayak  and  Dreger,  2014; 
Chiang  et  al.,  2014)  simple  but  wellcalibrated  ID  velocity  models  have  been  used  to 
compute  the  Green’s  function  and  invert  for  the  MT.  This  assumption  is  valid  at  relatively 
long  periods  where  the  wavelengths  being  inverted  are  much  larger  than  the  spatial 
dimension  of  the  sources,  and  the  source  time  history  is  short  compared  to  the  frequency 
passband  of  the  inversion.  However  the  limits  of  the  definitions  are  not  very  well  defined, 
nor  are  they  thoroughly  explored  in  the  published  literature. 

The  approach  to  compute  Green’s  functions,  especially  in  regions  of  low-seismicity 
where  high-resolution  velocity  models  are  not  available,  is  through  waveform  modeling  of 
regional  earthquakes  to  produce  calibrated  ID  velocity  models  (e.g.  Dreger  and 
Helmberger,  1990;  Bhattacharyya  et  al.,  1999;  Baise  et  al.,  2003;  Kim  et  al.,  2011). 
However,  the  ID  velocity  model  assumption  can  be  the  greatest  source  of  error  in  the  MT 
solution  and  has  never  been  thoroughly  investigated  for  the  source-type  discrimination 
application.  In  Ford  (2009b),  the  sensitivity  of  full  MT  solutions  for  the  2009  Memorial 
Day  DPRK  nuclear  explosion  to  ID  velocity  structure  was  examined.  In  that  study  880  ID 
velocity  models  for  the  region  from  the  North  Korean  test  site  to  regional  distance 
broadband  stations,  derived  from  the  Monte  Carlo  Markov  Chain  inversion  results 
(Pasyanos  et  al.,  2006),  were  tested  to  investigate  the  effect  on  the  recovery  of  the 
seismic  MT  source  type  (e.g.  Hudson  et  al.,  1989;  Ford  et  al.,  2010).  The  results 
indicated  that  while  there  was  a  strong  effect  on  the  source-type  and  the  ability  of  the 
Green’s  functions  for  the  individual  velocity  models  to  fit  the  data,  overall  the  velocity 
model  uncertainty  was  not  large  enough  to  affect  the  discrimination  of  the  seismic  event  as 
an  explosion.  The  spread  of  the  derived  solutions  in  source-type  space  (e.g.  Hudson  et  al., 
1989)  in  that  study  does  illustrate  however  that  the  effect  needs  to  be  examined  further, 
particularly  for  cases  where  precise  source-type  determination  is  important,  and  for  very 
sparse  monitoring  conditions  where  the  range  of  acceptable  solutions  in  terms  of  the 
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ability  to  fit  the  data  can  become  large  due  to  poor  station  coverage.  From  our  study  of 
explosions  in  Eastern  Kazakhstan  and  Lop  Nor  (Chapter  2;  Chiang  et  al.,  2014)  we  also 
observe  similar  variations  in  waveform  fits,  moment  and  MT  derived  source-types  as  a 
function  of  velocity  model.  Thus,  a  comprehensive  understanding  of  the  uncertainties 
associated  with  the  assumed  ID  velocity  structure  is  crucial  in  being  able  to  apply  the  MT 
based  discrimination  method  to  new  regions  of  interest  in  order  to  fully  assess 
populations  of  MT  results  for  explosions,  collapses,  earthquakes,  geothermal  and  volcanic 
events. 

There  are  some  studies  that  have  examined  errors  in  regional  MT  inversions  due  to 
biases  from  velocity  models.  Wei  et  al.  (2012)  shift  the  theoretical  seismograms  used  in 
their  ID  Green’s  functions  according  to  a  2D  correction  surface.  Comparison  with  and 
without  the  time  shift  shows  the  solution  improved  when  the  authors  incorporated  the 
correction.  This  is  similar  to  what  Ford  (2008)  found  in  which  when  time  shifts  were 
employed  the  spread  of  MT  solutions  for  different  ID  velocity  model  Green’s  functions 
in  source  type  space  was  reduced.  In  the  Berkeley  MT  code  (Dreger,  2003;  Minson  and 
Dreger,  2008)  a  time  shift  is  applied  for  optimal  alignment  of  the  data  with  the  applied 
Green’s  function  which  takes  into  account,  to  a  degree,  the  discrepancy  between  the 
velocity  model  and  the  actual  Earth  structure  (Pasyanos  et  al.,  1996).  Sileny  (2004)  has 
investigated  the  sensitivity  of  full  regional  MT  inversion  for  an  Mw4.8  earthquake  in  Italy 
to  perturbed  ID  velocity  models  and  found  that  the  DC  mechanism  is  well-resolved  from 
10  to  50  second  but  not  from  5  to  10  seconds.  Liu  et  al.  (2004)  obtain  MTs  for  three 
earthquakes  in  Southern  California  with  a  regional  3D  model  using  the  spectral  element 
method  and  compare  the  3D  solutions  with  those  obtained  from  body  wave  and  surface 
wave  inversions,  and  found  that  the  solutions  agree  well.  Covellone  and  Savage  (2012) 
calculated  3D  Green’s  functions  for  earthquakes  in  the  Middle  East  and  show  improved 
waveform  fits  and  a  reduction  of  non-double-couple  (non-DC)  components.  Panning  et 
al.  (2001)  included  near  source  3D  velocity  heterogeneity  and  found  its  effect  to  be 

minimal  in  the  bandwidth  of  interest,  although  small  isotropic  components  (~10%  of  the 

total  scalar  moment)  could  arise  due  to  the  unmodeled  near-source  3D  structure.  Although 
there  have  been  studies  on  the  effects  of  3D  structure  on  source  inversions  (e.g.  Panning 
et  al.,  2001),  to  our  knowledge,  there  has  not  been  much  work  investigating  the  effect  of 
3D  heterogeneity  along  the  full  source-receiver  path  on  the  full  regional  MT  for 
earthquakes  aside  from  the  Covellone  and  Savage  (2012)  work.  For  non-DC  sources, 
Tkalcic  et  al.  (2009)  determined  the  source  parameters  of  a  volcanic  earthquake  using  data 
from  20  to  50  seconds  and  path-specific  ID  velocity  models  to  account  for  3D 
heterogeneity.  This  work  offers  the  first  systematic  investigation  of  potential  bias  due  to 
unaccounted  3D  structure  on  regional  MT  inversion. 

As  a  result  of  the  unprecedented  coverage  of  the  NSF  EarthScope  transportable 
array  (TA),  a  number  of  high-resolution  regional  3D  Earth  models  have  been  published  in 
the  literature  (Bensen  et  al.,  2009;  Moschetti  et  al.,  2010;  Porritt  et  al.,  2014;  Shen  et  al., 
2013).  The  availability  of  these  finer  resolution  regional  tomography  models,  particularly 
the  ambient  noise  surface  wave  models  allow  us  to  generate  synthetic  waveforms  in  the 


Approved  for  public  release;  distribution  is  unlimited. 

41 


frequency  range  relevant  for  regional  waveform  modeling  of  small  to  moderately-large 
events.  Previous  work  by  Ford  et  al.  (2009a)  studied  a  variety  of  events  in  the  western  U.S. 
with  the  majority  of  events  near  or  at  the  Nevada  Test  Site  (NTS).  These  events  included 
nuclear  and  chemical  explosions,  earthquakes  and  mine  and  nuclear  cavity  collapses.  For 
these  events  we  study  the  impact  of  the  ID  velocity  model  assumption  on  moment  tensor 
results  be  means  of  a  series  of  synthetic  studies.  In  the  synthetic  studies  we  compute  3D 
velocity  model  GFs  using  the  anelastic  finite- difference  (FD)  code,  Seismic  Waves,  4th 
order  (SW4),  developed  at  the  Center  of  Applied  Scientific  Computing  at  Lawrence 
Livermore  National  Laboratory  (Sjogreen  and  Petersson,  2012;  Petersson  and  Sjogreen, 
2014).  We  use  the  Moschetti  et  al.  (2010)  tomographic  model  for  the  western  U.S.  and 
station  geometries  from  Ford  et  al.  (2009a).  We  then  invert  these  3D  synthetic  data  using 
GFs  for  a  ID  model  to  determine  the  seismic  MT  and  compare  the  results  to  the  known 
input  moment  MT.  The  effects  of  the  ID  velocity  model  assumption  can  then  be 
interpreted  via  the  difference  between  input  and  inverted  solutions. 

3.3.2  Velocity  Model  and  Computation  Setup 

The  synthetic  experiment  is  set  up  to  evaluate  the  effects  of  3D  heterogeneity  along 
the  source-receiver  path.  We  simulate  three-component  velocity  records  using  the 
Moschetti  et  al.  (2010)  3D  Earth  model  and  the  FD  code  SW4,  and  perform  the  standard 
MT  analysis  using  Green’s  functions  computed  from  a  well-calibrated  ID  Earth  model 
for  the  Western  U.S.  (Song  et  al.,  1996).  The  Moschetti  et  al.  (2010)  model  spans  the  entire 
western  U.S.  and  is  constructed  based  on  surface  wave  dispersion  measurements  from 
teleseismic  earthquakes  and  regional-scale  ambient  noise  measurements;  the  combination 
of  data  and  dense  station  coverage  better  constrains  the  shear  wave  velocity  structure  in  the 
crust  and  uppermost  mantle  than  previous,  earthquake  only  tomographic  models.  Moschetti 
et  al.  (2010)  invert  for  both  isotropic  and  anisotropic  models  but  in  this  study  we  only 
consider  the  isotropic  model.  The  3D  velocity  model  consists  of  a  low  velocity  sediment 
layer  and  three  crystalline  crustal  layers  above  the  Moho,  and  a  smoothly  varying  upper 
mantle  below  the  Moho.  The  compressional  wave  velocity  to  shear  wave  velocity  (Vp/Vs) 
ratio  and  density  structure  of  the  inverted  3D  model  are  based  on  empirical  relations 
between  wave  speed  and  density  (Brocher,  2005).  Focusing  on  California  and  Nevada 
(Figure  19  the  California  Coast  Ranges  have  slow  wave  speeds  in  the  upper  and  middle 
crust,  and  beneath  the  Central  Valley  of  California  the  low  wave  speeds  are  associated  with 
thick  sediments  in  the  San  Joaquin  Basin  to  the  south  and  the  Sacramento  Basin  in  the 
north.  In  the  Central  Valley  and  southernmost  end  of  Sierra  Nevada  the  Moho  varies 
between  20  to  40  km,  and  the  crust  thickens  near  the  eastern  edge  the  Sierra  Nevada. 
Moschetti  et  al.  (2010)  attribute  the  low  velocity  anomalies  of  the  lower  crust  beneath 
eastern  California,  the  Basin  and  Range  and  the  High  Lava  Plains  to  conductive  heating. 
The  Basin  and  Range  in  northern  Nevada  has  relatively  uniform  crustal  and  mantle 
structure,  as  the  relatively  long-period  surface  waves  used  to  construct  the  model  are  not 
sensitive  to  large  impedance  contrasts  from  very  thin  layers. 

For  our  FD  simulation  we  did  not  include  any  water  layers.  We  account  for 
attenuation,  Q,  using  a  function  of  local  Vs  [km/s]  (e.g.  Graves  et  al.,  2008;  Olsen  et  al., 
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2009),  such  that  Qs  =  50Vs  and  Qp  =  Qs.  The  FD  code  SW4  solves  the  seismic  wave 
equation  in  a  displacement  formulation  using  a  summation-by-parts  principle  and  it  is 
fourth-order  accurate  in  space  and  time  (Sjogreen  and  Petersson,  2012;  Petersson  and 
Sjogreen,  2014). 


Figure  19:  Background  color  represents  the  shear  wave  velocity  (Vs)  from  Moschetti  et  ah, 
(2010)  at  various  crustal  depths,  white  circle  is  the  event  location,  and  black  triangles  are 
the  station  locations. 


The  code  utilizes  a  Cartesian  mesh  to  solve  the  seismic  wave  equations  and  the  mesh 
can  be  built  using  the  PROJ.4  cartographic  projections  library  (Developed  by  G.  Evenden, 
https://trac.osgeo.org/proj/)  to  ensure  accurate  mapping  between  geographic  and  Cartesian 
coordinates.  In  contrast  to  a  staggered  grid  finite  difference  scheme  (e.g.  Larsen  and  Shultz, 
1995),  SW4  uses  a  node-centered  formulation  to  reduce  computation  memory.  The 
computation  domain  consists  of  a  free-surface  boundary  at  the  top  with  the  option  for 
irregular  topography  (Appelo  and  Petersson,  2009),  and  absorbing  super-grid  boundaries  in 
the  far  field  to  minimize  boundary  reflections  from  the  side  and  the  bottom  of  the 
computation  domain.  We  performed  the  calculations  on  the  Livermore  Computing  Center 
(LC)  8-core  Xeon  35-2670  Linux  cluster,  cab,  consisting  of  1,296  nodes,  each  with  16  cores 
per  node  and  32  GB  memory.  The  velocity  model  is  discretized  on  a  Cartesian  grid  with 
500  m  grid  spacing  and  a  minimum  Vs  of  1.1  km/sec.  SW4  internally  resamples  the  coarser 
3D  tomographic  model  (0.5  degree  grid  spacing)  on  to  the  finer  computation  grid  via  2- 
steps:  (1)  linear  interpolation  in  the  vertical  direction  and  (2)  Gaussian  averaging  in  the 
horizontal  directions.  The  model  domain  extends  approximately  665  km  x  950  km  x  70  km, 
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and  we  apply  a  Dirac  delta  forcing  function  at  the  source  to  simulate  a  seismic  source  in 
our  FD  calculations.  To  avoid  numerical  artifacts  due  to  grid  dispersion  we  consider  10  grid 
points  per  minimum  wavelength,  and  therefore  the  maximum  resolvable  frequency  is  0.2  Hz 
for  the  given  minimum  Vs.  Since  we  are  inverting  data  between  0.01  to  0.1  Hz  and  given 
the  magnitude  range  (M  <  5)  of  the  events  we  are  investigating  it  is  not  necessary  to 
consider  a  source  time  history  beyond  the  basic  Dirac  delta  function. 

Prior  to  performing  the  synthetic  sensitivity  tests,  we  performed  benchmark  testing  to 
validate  the  wave  propagation  simulated  using  SW4.  We  generated  and  compared  ID  velocity 
waveforms  using  the  two  different  methods:  (1)  the  frequency-wavenumber  (FK)  integration 
method  (Wang  and  Herrmann,  1980;  Herrmann  and  Wang,  1985;  Zhu  and  Rivera,  2002)  and 
(2)  the  numerical  FD  method  (SW4).  The  FK  code  is  included  as  part  of  a  software  package 
called  Computer  Programs  in  Seismology  (CPS3.30)  developed  by  Herrmann  (2013)  at  St. 
Louis  University.  We  use  a  simple  three-layer  model  with  no  attenuation  and  computed  the 
synthetics  at  100  km  distance  for  the  fundamental  Green’s  tensor  force  couples  (Figure  20 
Both  methods  agree  closely,  where  deviations  are  mostly  due  to  small  amplitude  differences 
in  the  Gxz ,  Gy  z  and  Gzz  components. 
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Figure  20:  Wavefonn  comparison  between  finite-difference  (SW4,  solid  black  line)  and 
frequency  wavenumber  integration  (FK,  dashed  line). 

3.3.3  Sensitivity  to  3D  Velocity  Structure 

We  set  up  the  synthetic  experiment  to  simulate  the  standard  processing  procedure 
used  in  regional  MT  analysis.  The  simulated  3D,  three-component  velocity  records  (north- 
south,  east-west  and  up-down)  from  SW4  are  treated  as  “data”  in  our  synthetic  test.  For  each 
station,  we  rotate  the  data  to  transverse,  radial  and  vertical,  and  add  random  Gaussian  white 
noise  to  the  data  by  scaling  the  amplitude  of  the  noise  to  10%  of  the  maximum  amplitude 
of  the  three  components.  ID  Green’s  functions  were  calculated  using  the  FK  method 
Herrmann  (2013).  The  3D  synthetic  data  and  ID  Green’s  functions  are  both  filtered  between 
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10  to  50  seconds.  Differences  between  the  input  and  recovered  MT  reveals  the  bias  due  to 
3D  structure.  A  direct  comparison  of  the  3D  and  ID  synthetics  for  a  normal  earthquake 
mechanism  shows  minor  differences  in  the  arrival  times  and  waveforms,  and  almost  no 
difference  in  the  20  to  50  second  period  range  which  is  a  commonly  used  passband  in 
regional  waveform  source  inversion  Pasyanos  et  al.  (1996).  This  is  due  to  the  large 
wavelengths  in  the  20  to  50  second  passband  and  the  smooth  nature  of  the  3D  tomographic 
model.  Nevertheless,  as  regions  of  monitoring  interest  are  calibrated  the  3D  models 
produced  are  likely  to  be  with  the  same  type  of  regional-scale  surface  wave  approach  used  to 
construct  the  3D  tomographic  model  examined  in  this  study.  The  following  section  describes 
the  results  for  the  three  types  of  sources  we  have  tested:  (1)  double-couple  (earthquake),  (2) 
isotropic  (explosion)  and  (3)  composite  source  mechanisms.  We  used  the  Little  Skull 
Mountain  mainshock  station  geometry  of  Ford  et  al.  (2009a)  for  all  of  the  synthetic  sources 
in  this  experiment  to  reflect  a  real  regional  monitoring  configuration. 

3.3.4  Double-Couple  Sources 

First  we  examine  the  ID  model  assumption  for  the  double-couple  (DC)  sources.  We 
tested  two  DC  mechanisms:  a  normal  earthquake  and  an  oblique  reverse  event,  both  at  4.5 
km  depth.  The  normal  faulting  mechanism  is  taken  from  the  Little  Skull  Mountain 
mainshock  MT  solution  in  Ford  et  al.  (2009a).  For  the  DC  synthetics,  both  deviatoric  and 
full  MT  inversions  can  recover  the  correct  source  mechanism  (Figure  21  and  the  variance 
reduction  (VR),  a  measure  of  the  goodness-of-fit  between  data  and  synthetics,  are  very 
similar  for  both  inversion  schemes.  Although  the  full  MT  results  have  up  to  16%  ISO 
component  in  the  full  MT  inversion  (Figure  21  the  change  in  VR  is  not  statistically 
significant  as  determined  by  an  F-test  (Menke,  1989;  Templeton  and  Dreger,  2006; 
Menke,  2012),  where  the  level  of  significance  is  only  57%.  Similarly,  the  Network 
Sensitivity  Solution  [NSS]  (Ford  et  al.,  2010),  a  mapping  of  the  goodness  of  fit  parameter  in 
the  complete  source-type  space,  shows  an  earthquake-like  surface  centered  on  the  origin 
(pure  DC).  The  CLVD  components  are  17%  and  40%  for  the  normal  and  oblique 
mechanisms,  respectively  (Figure  21  Interestingly  the  oblique  earthquake  was  found  to  have 
larger  variability  in  the  non-DC  components  compare  to  the  normal  earthquake  (Figure  22 
The  best  solutions  for  the  oblique  mechanism  are  more  widely  distributed  in  the  NSS.  The 
differences  in  the  NSS  distribution  and  non-DC  components  between  the  two  types  of  DC 
mechanisms  are  likely  the  result  of  different  radiation  patterns,  and  the  non-uniform 
sampling  of  the  radiation  from  the  sparse  station  coverage.  All  stations  except  one  are 
located  to  the  west  of  the  source  (Figure  19). 

Total  moment  from  the  full  MT  solution  is  quite  close  to  the  true  value,  estimates  of 
seismic  moment  are  all  within  10%  of  the  true  value.  The  source  depth  can  be  determined 
from  the  goodness  of  fit  curve  (Figure  22  however  there  is  some  mechanism  sensitivity.  The 
source-type  plots  (Figure  22  illustrate  the  tradeoff  between  source  depth  and  mechanism, 
shwoing  that  as  source  depth  deviates  from  the  true  depth  the  ISO  component  also  increases. 
But  for  both  cases  the  MT  solutions  are  stable  with  respect  to  source  depth  and  the 
solutions  with  the  best  waveform  fits  (VR  >  60%)  are  between  5  and  7  km.  The  best  fitting 
solutions  yield  the  correct  mechanism  and  depths  close  to  the  correct  source  depth  of  4.5 
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km.  In  summary  for  the  two  DC  cases,  we  can  recover  the  correct  fault  plane  solution  using 
only  regional  waveforms  and  it  appears  that  the  ID  model  assumption  is  adequate,  but  one 
must  be  careful  when  trying  to  interpret  the  CLVD  component.  Apparently  the  effect  of  3D 
heterogeneity  is  negligible  at  these  periods  for  recovering  the  correct  DC  mechanism  using  a 
ID  model,  but  3D  heterogeneity  can  increase  contributions  from  the  CLVD.  However,  the 
NSSs  show  a  pure  DC  mechanism  fits  the  data  equally  well  and  can  be  used  to  quantify  the 
level  of  variability  in  the  non-DC  components  (Fig.  11a). 
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Figure  21 :  The  simulated  3D  data  is  in  black,  the  deviatoric  solution  is  in  green  and  the 
full  solution  is  in  brown.  All  waveforms  are  in  displacement  (cm). 
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(b)  Normal  Oblique 


Figure  22:  (a)  Network  Sensitivity  Solutions  color-coded  by  variance  reduction  (VR) 
presented  on  the  Tape  and  Tape  (2012a)  and  Tape  and  Tape  (2012b)  Lune.  The  white  circle 
represents  the  location  of  the  best  full  moment  tensor  solution  in  source-type  space,  (b)  VR 
with  respect  to  source  depth  for  two  different  DC  mechanisms  and  color-coded  by  source 
depth.  The  filled  circles  are  full  moment  tensor  solutions  and  the  filled  squares  are  deviatoric 
moment  tensor  solutions.  The  full  moment  tensor  solutions  are  also  plotted  in  source-type 
space  and  color-coded  by  source  depth  as  well. 

3.3.5  Isotropic  Sources 

A  direct  comparison  between  ID  and  3D  synthetic  data  for  the  explosion  case  shows 
more  variation  in  waveforms  compared  to  the  earthquake  cases  (Figure  23  which  is  likely 
due  to  the  shallow  explosion  source  depth  that  was  considered  (500  m).  In  addition,  the 
averaged  ID  model  is  less  representative  of  the  shallow  structure  in  the  3D  tomographic 
model.  Since  the  isotropic  source  radiation  pattern  has  no  azimuthal  dependence,  variations 
in  waveforms  due  to  the  3D  structure  are  more  prominent;  whereas  an  earthquake  source 
imposes  an  order  of  magnitude  variation  in  the  radiation  pattern  as  a  function  of  azimuth.  At 
long  periods  the  source  variation  is  large  compare  to  the  effects  of  3D  structure  in  this  long- 
period  passband.  The  recoverability  of  the  explosion  source  depth  is  similar  to  what  was 
found  in  the  previous  ID  synthetic  studies  (Fig.  11)  and  real  explosions  (Chiang  et  ah, 
2016)  in  that  the  source  depth  recoverability  primarily  depends  on  the  frequency  band  used 
in  the  MT  inversion.  Because  the  source  depth  is  much  less  than  the  wavelengths 
considered,  we  cannot  resolve  the  source  depth  in  detail  and  can  only  constrain  it  to  be 
shallower  than  1  km.  The  MT  analysis  shows  the  deviatoric  MT  solutions  are  stable  at 
shallow  depths  (Figure  24c)  but  the  isotropic  component  is  not  well  constrained  in  the  full 
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MT  inversion  since  the  isotropic  component  is  very  sensitive  to  the  time  shifts  in  the 
inversion.  The  full  MT  solutions  alternate  between  an  explosion  and  a  CLVD  mechanism 
with  the  major  vector  dipole  in  compression  depending  on  the  time  shifts  due  to  the  fact 
that  the  long-period  surface  wave  radiation  of  the  isotropic  source  and  a  CLVD  with  a  major 
vector  dipole  in  compression  are  the  same.  The  tradeoff  between  these  two  source  types  is 
seen  in  the  NSS  for  this  test  (Figure  24c)  where  a  pure  explosion  mechanism  obtained  from 
a  grid-search  method  has  a  VR  of  75%  and  the  full  MT  has  a  VR  of  76%. 

The  alternation  between  the  two  mechanisms  is  the  manifestation  of  the  classic 
tradeoff  for  an  explosion  source  (e.g.  Ford  et  ah,  2009b).  For  an  isotropic  source,  the  level 
of  false  non-ISO  components  in  the  moment  tensor  results  from  unmodeled  3D  structures 
and  the  theoretical  ISO-CLVD  tradeoff,  and  thus  it  is  difficult  to  tease  out  how  much  of  the 
CLVD  component  is  due  to  unmodeled  3D  structure.  Combining  waveform  and  first  motion 
polarities  can  eliminate  the  tradeoff  to  get  the  best  MT,  and  we  obtain  a  solution  located 
close  to  the  theoretical  explosion  (Figure  24c).  The  best  solution  consists  of  a  predominantly 
explosive  mechanism  with  60%  ISO,  28%  CLVD  and  12%  DC.  The  MT  still  has  a  relatively 
high  CLVD  component  but  the  level  of  false-DC  is  low.  This  is  an  important  result  for 
discrimination  because  it  shows  that  unmodeled  3D  structure  at  the  typical  passband  used 
in  regional  distance  nuclear  monitoring  does  not  introduce  false-DC  in  the  inversion  results. 
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Figure  23:  Dashed  lines  are  waveforms  calculated  using  a  ID  Earth  model  and  solid  lines 
are  waveforms  calculated  using  a  3D  Earth  model.  The  waveforms  are  in  displacement 
(cm)  and  band-pass  filtered  between  10-50  seconds. 
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a.  Moment  tensor  inversion  results 
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b.  VR  as  a  function  of  source  depth 
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Figure  24:  (a)  Explosion  moment  tensor  solutions  and  waveform  fits  from  10  to  50  seconds. 
The  simulated  3D  data  is  in  black,  the  deviatoric  solution  is  in  green  and  the  full  solution  is  in 
brown.  All  waveforms  are  in  displacement  (cm),  (b)  Variance  reduction  (VR)  as  a  function 
of  source  depth.  Circles  are  full  moment  tensor  solutions  and  the  squares  are  deviatoric 
moment  tensor  solutions,  (c)  Network  Sensitivity  Solutions  using  only  waveforms,  only  first 
motions  or  combining  waveforms  and  first  motions.  The  white  circle  represents  the  location 
of  the  best  full  moment  tensor  solution  in  source-type  space. 
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3.3.6  Composite  Sources 

Explosion  processes  are  often  not  purely  isotropic,  and  they  typically  exhibit  evidence  of 
complex  seismic  radiation,  particularly  from  the  presence  of  long-period  Love  wave  energy. 
Large  non-ISO  components  have  been  observed  in  real  explosion  data,  sometimes  referred  to 
as  tectonic  release  (e.g.  Masse,  1981;  Ekstrom  and  Richards,  1994),  and  this  component  can 
become  large  enough  to  reverse  the  Rayleigh  wave  polarities  (e.g.  Walter  and  Rodgers,  1999; 
Chiang  et  ah,  2014)  which  can  complicate  discrimination,  as  well  as  generate  long-period  (20 
to  50  s)  Love  waves  that  can  have  the  largest  amplitudes  of  all  three-components  at  a  given 
station  (e.g.  Dreger  and  Woods,  2002;  Lord  et  al.,  2009a).  In  addition,  shear  failure  is  often 
seen  in  real  explosions  (Toksoz,  et  al.,  1965;  Toksoz  and  Kehrer,  1972;  Burger  et  al.,  1986; 
Day  et  al.,  1987).  Therefore  it  is  necessary  to  examine  whether  the  MT  inversion  method 
breaks  down  and  can  become  biased  when  both  the  effects  of  complex  source  processes  and 
unaccounted  for  3D  effects  are  included,  and  whether  the  ability  to  successfully  identify  a 
given  seismic  event  as  an  explosion  can  be  compromised. 

The  composite  source  considered  here  consists  of  a  pure  explosion  and  reverse 
mechanism  (strike=  135;  dip=  40;  rake=  105)  both  at  500  m  depth.  We  apply  the  same  data 
processing  procedure  and  station  configuration  used  for  the  DC  and  ISO  cases,  and  added 
random  Gaussian  white  noise.  For  a  composite  source  that  is  predominantly  explosive,  the 
combination  of  3D  effects  and  source  tradeoff  appears  to  have  an  impact  on  the  inverted 
solution.  The  level  of  CLVD  can  be  as  high  as  60%  (for  the  case  with  40%  DC).  The  MT 
solutions  and  waveforms  fits  for  a  composite  mechanism  with  30%  DC  is  shown  in  Figure 
25a.  The  full  MT  solution  is  45%  CLVD  suggesting  that  the  solution  is  still  dominated 
by  the  theoretical  trade-off  for  explosive  mechanism.  Figure  25b  shows  that  the  solutions  for 
composite  sources  with  varying  percentages  of  DC  are  still  dominated  mostly  by  the  explosion- 
CLVD  trade-off  and  the  waveform  NSS  distributions  remain  more  explosive-like  even  as  the 
contribution  from  the  DC  increases.  Although  the  recovered  DC  components  for  the  MT 
inversions  are  close  to  the  input  value  (30%  DC),  the  inversions  do  not  recover  the  correct  DC 
mechanism.  Instead  of  a  reverse  mechanism  the  full  MT  inversion  for  a  composite  source 
with  30%  DC  recovers  a  normal  mechanism  (strike=  169;  dip=  72;  rake=  —104). 

The  inversion  recovers  a  shallow  dipping  reverse  component  only  when  the  DC 
contribution  reaches  50%  (Figure  25b).  The  orientation  of  the  DC  cannot  be  resolved  for  a 
predominately  explosive  composite  source.  The  results  of  this  test  and  the  previous  tests 
indicate  that  there  can  be  substantial  variability  in  the  minor  components  of  a  recovered 
moment  tensor.  These  minor  components  appear  to  be  the  most  strongly  affected  from  the 
unaccounted  for  3D  velocity  model.  It  is  good  news,  however  that  the  primary  moment 
tensor  components,  namely  a  DC  for  tectonic  earthquakes,  and  the  isotropic  component  for 
explosions  are  not  strongly  affected  and  therefore  the  discrimination  of  the  two  types  of 
mechanism  remains  robust. 
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(b)NSS 


(a)  Composite  source  with  30%  doube-couple  (DC) 

Tangential  Radial  Vertical 


DC  =  28% 


Figure  25:  (a)  Deviatoric  (green)  and  full  (brown)  moment  tensor  solutions  and  waveform 
fits.  The  simulated  data  for  an  explosion  plus  reverse  fault  mechanism  are  in  black.  The 
example  shown  here  is  a  composite  source  with  30%  DC.  (b)  The  best-fitting  full  moment 
tensor  solution  and  NSS  in  source-type  space  for  four  composite  sources.  All  solutions 
have  the  same  mechanism  (explosion  +  reverse  fault)  but  different  percentage  of  DC.  The 
contours  represent  moment  tensor  solutions  with  normalized  VR>  98%. 


3.3.7  Path  Calibrations 

One  solution  for  trying  to  reduce  possible  MT  solution  bias  due  to  unaccounted  for 
3D  wave  propagation  effects  is  to  develop  source-receiver  specific  ID  velocity  models. 
Given  earthquake  data  in  the  source  region  and  independent  constraint  on  the  focal 
mechanism,  it  would  be  possible  to  perform  broadband  waveform  modeling  either  by 
forward  modeling,  genetic  algorithm,  or  full  grid-search  approaches  (e.g.  Dreger  and 
Helmberger,  1990;  Bhattacharyya  et  ah,  1999;  Baise  et  al.,  2003;  Kim  et  al.,  2011)  to 
improve  the  source-receiver  velocity.  On  the  other  hand  if  there  is  independent  information 
from  tomography  and  receiver  functions  the  information  can  be  used  to  compute  path 
specific  ID  models  to  improve  seismic  MT  recovery  in  regions  of  large  lateral  heterogeneity 
(Tkalcic  et  al.,  2009).  Here  we  use  the  3D  synthetic  data  to  evaluate  the  procedure  of  using 
source-receiver  path  specific  ID  velocity  structure  extracted  from  tomographic  models. 

We  use  the  averaged  elastic  properties  and  attenuation  from  the  3D  tomographic  model  for 
each  source-receiver  pair,  with  a  total  of  seven  path-specific  ID  velocity  models  (Figure 
26a)  and  compute  the  Green’s  functions  for  each  of  these  models.  The  average  ID  models 
have  slightly  shallower  sediment  layer  and  mantle  velocities  compared  to  the  Song  et  al. 
(1996)  model.  But  the  Moho  depths  are  very  similar  (within  5  km)  for  all  the  models. 
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Our  synthetic  study  shows  path-specific  calibration  does  help  in  cases  where  one  simple  layered 
model  cannot  adequately  characterize  the  underlying  Earth  structure  between  different  source- 
receiver  pairs,  it  reduces  the  non-DC  components  in  earthquakes  in  the  time-domain  wavefonn 
inversion  and  increases  the  waveform  fits.  The  MTs  are  88%  and  69%  DC  for  the  normal  and 
oblique  mechanisms,  respectively.  The  biggest  difference  in  the  NSS  using  a  universal  ID 
model  for  all  stations  and  using  path-specific  ID  models  for  the  pure  DC  sources  is  that  for 
the  path-specific  models  the  goodness  of  fit  distribution  is  centered  on  the  pure  DC 
solution  and  is  not  shifted  along  the  deviatoric  or  isotropic  axes.  Comparing  Figure  26b 
to  Figure  22  we  see  that  the  region  of  best  fitting  solutions  (normalized  VR>  98%)  has  a 
smaller  distribution  for  the  oblique  mechanism  whereas  for  the  normal  mechanism  the 
area  of  the  contoured  region  remains  similar  to  the  results  using  a  single  velocity  model. 
For  the  explosion  case,  the  tradeoff  between  explosion  and  CLVD  persists  although  the  best 
solution  becomes  more  isotropic  compared  to  the  result  without  path  calibration.  For  the 
composite  source  with  30%  DC,  the  best  MT  is  explosive  and  the  CLVD  has  reduced  to 
only  12%.  However  we  still  cannot  recover  the  correct  DC  mechanism,  which  is  surprising. 
The  recovered  fault  plane  solution  is  a  shallow  dipping  normal  mechanism.  In  this  study, 
we  found  that  the  DC  sources  benefited  more  from  the  path-specific  approach. 


Figure  26:  (a)  Seven  ID  models  derived  from  the  3D  model  by  Moschetti  et  al.,  (2010). 
(b)  Full  moment  tensor  solutions  and  NSS  result  using  path-specific  ID  models  for  each 
source-receiver  pair.  The  contours  are  the  best-fitting  solutions  with  a  normalized 
VR>  98%. 
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3.3.8  Random  Velocity  Perturbations 

Since  the  3D  tomographic  models  are  smooth  representations  of  Earth  structure  we 
wanted  to  examine  a  more  severe  case  in  which  we  apply  random  velocity  perturbations 
in  the  model  to  increase  the  effects  of  scattering  and  multipath  in  the  3D  synthetic  data. 
A  realistic  description  of  the  inhomogeneous  Earth  can  be  described  by  the  von  Karman 
correlation  function  (Goff  and  Jordan,  1988): 


(26) 


Ku  is  the  modified  Bessel  function  of  the  second  kind  order  u,  where  0  <  u  <  1 .  u  is  the 
Hurst  number  which  describes  the  roughness  of  the  medium,  T  is  the  Gamma  function,  x,  y, 
and  z  are  vectors  in  the  3D  random  field,  and  ax,  ay ,  az  are  the  correlation  distances.  The 
von  Karman  correlation  function  is  self- variant  and  rich  in  short  wavelength  components, 
making  it  a  common  approach  to  model  geologic  structures  with  a  fractal  nature  over  a  wide 
scale  (Wu,  1982).  The  Hurst  number,  which  controls  the  correlation  decay,  is  typically  less 
than  0.5  from  modeling  of  seismic  reflection  data  (e.g.  Holliger  and  Levander,  1992;  Nielsen 
and  Thybo,  2006).  Correlation  distances  on  the  order  of  a  few  kilometers  to  tens  of 
kilometers  have  minimal  effect  on  the  10  to  50  second  waves,  therefore  very  large  values  are 
necessary  to  observe  changes  in  the  synthetic  waveforms  relative  to  the  smooth  3D  model 
in  the  period  range  of  interest.  We  use  a  Hurst  number  of  0.3  and  horizontal  and  vertical 
correlation  lengths  of  250  km  and  50  km,  respectively,  to  perturb  the  crustal  velocities.  We 
apply  velocity  fluctuations  up  to  20%  (Figure  27)  and  obtain  a  rough  model  with  crustal 
velocities  following  a  von  Karman  distribution.  Figure  28a  shows  the  perturbed  model  at  5 
km  depth.  The  velocity  contrast  between  the  California  coast  and  the  Basin  and  Range  is 
preserved  but  with  rough  edges.  Figure  28b  shows  a  cross  section  of  crustal  velocities  from 

the  coast  to  the  Nevada  Test  Site,  near  the  location  of  the  synthetic  source  (  —  116 

degree  longitude).  The  layers  are  well  defined  but  with  small-scale  heterogeneities  scattered 
throughout  the  crust. 

We  generate  3D  synthetics  using  the  perturbed  model  for  a  normal  earthquake  and 
an  explosion,  and  compare  the  MTs  to  the  solutions  from  the  smooth  model.  A  comparison 
between  normal  earthquake  synthetics  generated  from  the  ID  model  and  perturbed  3D 
model  low-pass  filtered  at  0.15  Hz  shows  the  complexities  due  to  scattering  in  the  3D 
synthetics  (Figure  29)  The  scattering  (ringing)  on  the  tangential  component  in  the  ID 
synthetics  is  due  to  the  low  velocity  layer  in  the  ID  model  (Fig.  26a).  When  these  data 
were  inverted  from  10  to  50  seconds  a  normal  mechanism  was  obtained  but  the  goodness 
of  fit  was  lower  than  the  typical  level  (VR>  60%)  considered  indicating  a  well-constrained 
solution.  Figure  30a  shows  the  deviatoric  MT  is  83%  DC,  which  is  10%  lower  than  the 
solution  using  the  smooth  model,  whereas  the  full  MT  remains  similar  to  the  previous 
result  (Figure  21  with  68%  DC.  But  the  NSS  (Fig.  30b)  shows  the  population  of  best- fitting 
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solutions  (VR>46%)  deviates  further  away  from  a  pure  DC  source  and  are  more  broadly 

distributed  compare  to  Figure  22  The  increase  in  the  area  of  the  best- fitting  MTs  can  be 
interpreted  as  increasing  bias  in  the  solution  due  to  unmodeled  structure. 

In  contrast  to  results  for  the  smooth  3D  model,  explosion  deviatoric  and  full  MTs 
(Fig.  31a)  show  false-DC  components  increased  to  almost  50%.  VR  is  57%  for  both 
deviatoric  and  full  MTs,  a  value  close  to  the  level  of  fit  for  a  well-constrained  solution,  and 
a  pure  explosion  mechanism  obtained  using  a  grid-search  method  has  a  VR  of  52%. 
Although  the  waveform  NSS  (Fig.  31b)  still  shows  an  explosion-like  distribution  the  best- 
fitting  solutions  encompass  almost  half  of  the  area  in  the  source-type  space,  including  a 
pure  DC  source.  The  first  motion  constrained  NSS  eliminated  the  tradeoff  and  pure  DC 
source  types  but  compared  to  Figure  24c  the  uncertainties  estimated  from  the  constrained 
NSS  are  slightly  larger.  The  first  motion-constrained  full  MT  solutions  for  the  smooth  and 
perturbed  3D  model  are  both  60%  ISO,  but  the  solution  for  the  perturbed  model  has  a 
slightly  greater  DC  of  31%.  In  this  study  we  show  that  the  MT  uncertainty  increases 
when  the  3D  model  becomes  more  heterogeneous  and  the  false-DC  component  increases 
when  an  explosion  source  is  considered. 


Velocity  Perturbations 


Figure  27:  Velocity  fluctuations  following  a  von  Karman  distribution  with  a  Hurst 
number  of  0.3,  horizontal  and  vertical  correlation  distances  of  250  km  and  50  km, 
respectively. 
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3.3.8  Conclusion 


The  results  of  the  synthetic  3D  model  sensitivity  tests  show  that  using  ID  velocity 
models  to  compute  Green’s  functions  for  use  in  long-period  (as  short  as  10  seconds)  MT 
inversions  is  reasonable,  and  that  at  least  synthetically  we  do  not  find  significant  bias  in 
solutions,  nor  issues  in  being  able  to  discriminate  different  types  of  sources.  For  the  DC 
case,  we  can  recover  the  correct  mechanism  and  source  depth  using  a  universal  ID  velocity 
model  for  all  stations.  The  ISO  component  in  the  full  MT  is  less  than  20%  for  the  two  DC 
mechanisms  explored.  But  for  the  oblique  mechanism  the  inversion  recovers  a  solution  with 
a  highCLVD  component,  resulting  in  a  solution  that  is  50%  DC.  The  observation  of  false 
non-DC  components  up  to  50%  for  the  earthquake  source  shows  that  although  the  correct 
fault  plane  solutions  for  the  DC  can  be  recovered  at  relatively  long  periods,  the  CLVD 
component  can  be  large  in  full  MT  inversions,  however  the  NSS  remains  earthquake-like 
focused  at  the  origin  of  the  source-type  plot.  For  the  explosion  cases,  the  level  of  false  non- 
ISO  is  dominated  by  the  1SO-CLVD  tradeoff,  although  the  recovered  MT  consists  of  a 
CLVD  solution  the  DC  component  is  less  than  10%,  and  the  NSS  exhibits  the  typical 
explosion-like  signature  in  source-type  space  (e.g.  Ford  et  al.,  2010).  Including  additional 
data  from  first  motion  polarities  can  eliminate  the  trade-off  and  the  correct  solution  can  be 
recovered.  The  constrained  MT  solution  is  predominantly  explosive  and  has  a  very  low  DC 
of  9%.  Similarly,  results  for  the  explosive  composite  source  show  the  tradeoff  affects  the 
recoverability  of  the  MTs.  The  fault  plane  orientation  of  the  DC  component  in  a  composite 
source  cannot  be  constrained  well,  except  for  the  case  with  50%  DC  and  50%  ISO.  Results 
from  synthetic  tests  using  a  more  heterogeneous  3D  model  with  crustal  velocity  fluctuations 
following  the  von  Karman  distribution  show  the  full  MT  inversion  can  recover  the  correct 
DC  mechanism.  The  percentage  of  non-DC  remains  the  same  compared  to  the  result  for 
a  smooth  model  but  the  solution  has  higher  uncertainties  as  indicated  by  lower  waveform 
fits  and  larger  variability  in  the  NSS.  When  an  explosion  source  is  considered  the  false-DC 
component  increases  to  almost  50%  in  the  inversion.  The  NSS  exhibits  an  explosion-like 
signature  in  the  source-type  space  but  with  large  variability  that  extends  to  a  pure  DC  source. 
But  considering  the  best  solution  with  first  motion  constraints  shows  a  predominantly 
explosive  source  with  a  slightly  higher  DC  of  28%  compare  to  the  results  from  a  smooth 
model.  The  false  non-DC  components  in  earthquakes  and  the  false-DC  components  in 
explosions  arising  from  unaccounted  for  3D  path  effects  can  be  reduced  by  using  path- 
specific  ID  Green’s  functions.  The  Green’s  functions  are  computed  by  taking  the  average 
ID  model  for  each  source-receiver  path  from  the  3D  tomographic  model.  Path-calibration 
also  reduces  the  large  uncertainty  in  the  NSS  for  the  oblique  mechanism. 
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Depth  (km) 


Figure  28:  (a)  Shear  wave  velocities  (Vs)  at  5  km  depth.  Black  line  denotes  the  location 
of  the  cross  section,  (b)  Cross  section  of  the  crustal  velocities  at  depth. 
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ID  3D 


Figure  29:  ID  and  3D  synthetic  waveforms  for  a  normal  mechanism  (strike=202; 
rake=— 100;  dip=36)  at  SAO.  3D  synthetics  are  computed  from  the  perturbed  model. 


Figure  30:  (a)  Normal  earthquake  moment  tensor  solutions  and  waveform  fits  at  10  to 
50  seconds.  The  simulated  3D  data  using  the  perturbed  model  is  in  black,  the  deviatoric 
solution  is  in  green  and  the  full  solution  is  in  brown.  All  waveforms  are  in  displacement 
(cm),  (b)  NSS  using  only  waveforms  and  combined  waveforms  and  first  motions. 
Solutions  are  color-coded  by  VR.  The  white  circle  represents  the  location  of  the  best  full 
moment  tensor  solution  in  source-type  space. 
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b.  Network  Sensitivity  Solutions  (NSS) 
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Figure  31:  (a)  Explosion  earthquake  moment  tensor  solutions  and  waveform  fits  at  10  to  50 
seconds.  The  simulated  3D  data  using  the  perturbed  model  is  in  black,  the  deviatoric 
solution  is  in  green  and  the  full  solution  is  in  brown.  All  waveforms  are  in  displacement 
(cm),  (b)  NSS  using  only  waveforms  and  combined  waveforms  and  first  motions.  Solutions 
are  color-coded  by  VR.  The  white  circle  represents  the  location  of  the  best  full  moment 
tensor  solution  in  source-type  space. 
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3.4  Application  of  3D  Green’s  Functions  on  Recovery  of  NTS  Explosion 
Moment  Tensors 

3.4.1  Introduction 

Obtaining  reliable  source  mechanisms  of  small  magnitude  seismic  events  at  regional 
distance  can  be  challenging  due  to  poor  station  coverage,  high  noise  levels  in  the  long 
period  waves,  or  both.  Commonly  long  period  waveforms  are  utilized  to  minimize  errors 
due  to  unaccounted  Earth  structure,  but  with  small  magnitude  events  the  signal-to-noise 
ratios  (SNR)  of  long  period  waves  often  degrade  rapidly  with  increasing  distance  from  the 
source.  The  need  for  regional  distance  monitoring  of  low  yield  explosions  therefore 
necessitates  the  use  of  shorter  period  information  that  may  have  higher  SNR,  but  which 
also  requires  Green’s  functions  (GFs)  that  account  for  the  more  complex  short  period  wave 
propagation.  Recent  advancements  in  waveform  modeling  and  high-performance  computing 
techniques  have  made  the  use  of  source  inversions  at  relatively  short  periods  and 
structurally  complex  regions  an  attractive  and  viable  option,  and  has  been  shown  in  several 
studies  (Ramos-Martinez  and  McMechan,  2001;  Liu  et  ah,  2004;  Zhao  et  al.,  2006;  Walter 
et  al.,  2010;  Covellone  and  Savage,  2012)  to  be  effective.  However,  there  has  not  been  a 
systematic  analysis  of  the  impact  of  3D  velocity  structure  on  the  estimation  of  seismic 
source  parameters  for  the  purpose  of  explosion  monitoring. 

The  efficacy  of  GFs  depends  on  the  accuracy  of  the  3D  model.  Using  a  3D 
reference  model  to  account  for  complex  wave  propagation  in  and  near  the  Los  Angeles 
Basin,  Liu  et  al.  (2004)  and  Lee  et  al.  (2011)  computed  source  mechanisms  of  small  to 
moderate-sized  earthquakes  in  Southern  California  using  the  spectral  element  method  and 
finite-difference  method,  respectively,  and  found  good  agreement  between  the  ID  and  3D 
solutions.  Hingee  et  al.  (2011)  observe  improvements  in  3D  synthetics  over  ID  synthetics 
for  some  regions  in  Australia.  The  reduction  in  waveform  misfit  is  limited  to  certain 
regions  suggesting  the  3D  model  for  the  Australian  region  needs  further  improvement. 
Covellone  and  Savage  (2012)  computed  deviatoric  and  full  moment  tensors  (MTs)  for 

195  Mw~5.5  earthquakes  in  the  Middle  East  and  found  little  difference  among  the 

solutions  due  to  the  fact  that  a  large  number  of  observations  are  available  to  provide 
good  constraints  in  the  full  MT  inversion.  They  compared  the  deviatoric  ID  and  3D  MTs 
and  found  that  there  is  an  improvement  in  waveform  fits  and  a  reduction  in  the  non- double¬ 
couple  (non-DC)  components  when  using  3D  GFs.  Although  it  is  important  to  note  that  in 
their  comparison  study,  the  3D  synthetics  were  evaluated  against  ID  synthetics  computed 
from  the  preliminary  reference  Earth  model  (PREM;  Dziewonski  and  Anderson,  1981). 

The  focus  of  most  studies  has  been  on  earthquakes,  the  application  of  3D  source 
inversions  on  non-DC  events  have  not  been  thoroughly  investigated.  Walter  et  al.  (2010) 
computed  full  MT  solutions  for  icequakes  near  the  base  of  the  Swiss  Alpine  Glacier  and 
determined  seismic  events  are  results  of  near-horizontal  tensile  cracks.  The  authors 
concluded  the  ID  models  provided  robust  estimates  of  the  source  parameters,  although  the 
more  complex  3D  model  improved  the  fits  between  data  and  synthetics,  particularly  for 
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reflections  from  the  basal  moraine  and  rock.  There  has  not  been  a  systematic  study  of  the 
effects  of  3D  GF  on  the  recovery  of  seismic  MT  of  explosions  and  the  effect  on 
discrimination  capability. 

The  purpose  of  this  study  is  to  develop  a  platform  to  perform  routine  3D  MT 
inversions  and  to  compare  source-type  discrimination  results  for  explosions  and 
earthquakes  using  both  ID  and  3D  GFs.  The  same  techniques  (e.g.  Ford  et  ah,  2010) 
developed  to  evaluate  the  quality  of  ID  models  can  also  be  applied  to  3D  models.  Ford  et 
al.  (2009a)  obtained  full  MT  solutions  using  a  well-calibrated  ID  model  (Song  et  al.,  1996) 
for  explosions,  earthquakes  and  mine  collapses  in  the  western  United  States.  In  this  study 
we  compare  ID  and  3D  full  MT  solutions  at  two  frequency  bands:  20  to  50  seconds  and  8 
to  20  seconds.  The  high  frequency  cutoff  is  limited  by  the  resolution  of  the  Earth  model. 
We  use  the  3D  surface  wave  tomography  model  by  Moschetti  et  al.  (2010)  determined 
from  surface  wave  dispersion  measurements  at  periods  as  short  as  6  seconds.  Of  the  32 
explosions,  earthquakes  and  collapses  previously  studied  by  Ford  et  al.  (2009a),  we  select 
six  explosions  and  earthquakes  in  the  vicinity  of  the  former  Nevada  Test  Site  (NTS)  for 
source  inversions.  The  events  are  listed  in  Table  3  where  the  earthquake  event  information 
comes  from  the  Northern  California  Seismic  Network  (NCSN)  and  the  explosion  event 
information  comes  from  Springer  et  al.  (2002).  We  apply  the  source-receiver  reciprocity 
principle  to  compute  3D  GFs  using  the  finite-difference  (FD)  method  (Eisner  and  Clayton, 
2001;  Graves  and  Wald,  2001).  The  advantage  of  seismic  reciprocity  is  a  drastic  decrease  in 
computation  cost,  especially  when  the  number  of  sources  outweighs  the  number  of 
receivers,  because  for  a  set  of  three  FD  simulations  in  which  unit  forces  are  applied  at  the 
receiver  location  in  the  north,  east  and  vertical  directions  the  complete  MT  GF  can  be 
determined  at  all  interior  points  in  the  model.  Thus  in  the  monitoring  of  a  region  it  is 
possible  to  develop  GF  over  the  desired  depth  range  for  a  calibrated  3D  model  obtained 
independently  from  tomography  or  other  methods. 


Table  3:  Event  Information 


Name 

Date 

Latitude 

Longitude 

Depth  (km) 

Magnitude 

Amargosa 

1996/09/05, 

08:16:56.09 

36.6827 

-116.3378 

5000 

3.38  (Md) 

Little  Skull  Mt. 

2002/06/14, 

12:40:45.82 

36.6438 

-116.3448 

8750 

4.32  (Md) 

METROPOLIS 

1990/03/10, 

16:00:00.08 

37.112 

-116.056 

469 

4.94  (Md) 

COSO 

1990/03/10, 

16:00:00.08 

37.104 

-116.075 

417 

4.50  (Md) 

HOYA 

1991/09/14, 

19:00:00.08 

37.226 

-116.429 

658 

5.40  (Md) 

JUNCTION 

1992/03/26, 

16:30:00.00 

37.272 

-116.361 

622 

4.82  (Md) 
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3.4.2  Method 


We  perform  the  MT  inversion  in  the  time-domain  using  three-component 
displacement  seismograms  for  receivers  located  in  the  western  United  States  (Figure 
32).  The  data  are  instrument  corrected,  rotated  to  the  radial,  transverse  and  vertical 
components  and  bandpass  filtered  between  8  to  50  seconds  period  using  an  acausal  4-pole 
Butterworth  filter.  For  the  3D  MT  inversions,  an  additional  step  is  taken  by  convolving 
the  data  with  the  Gaussian  source  time  function  used  to  compute  the  3D  GFs  via  the  FD 
method.  This  is  done  since  it  is  more  stable  to  convolve  the  finite-difference  source  time 
function  with  the  data,  essentially  applying  a  filter,  than  it  is  to  deconvolve  the  source 
time  function  from  the  computed  Green’s  functions.  We  use  data  recorded  at  stations 
from  the  Berkeley  Digital  Seismic  Network  (BSDN),  Southern  California  Seismic 
Network  (SCSN)  and  the  Lawrence  Livermore  Seismic  Network.  The  United  States 
Geological  Survey,  Berkeley  Seismological  Laboratory  and  California  Institute  of 
Technology  jointly  manage  the  networks  in  California.  BDSN  and  SCSN  data  and 
instrument  responses  can  be  obtained  from  the  Incorporated  Research  Institutions  for 
Seismology  Data  Management  Center  (IRIS  DMC,  http://ds.iris.edu/ds/nodes/dmc/) 
and  the  Northern  California  Earthquake  Data  Center(http://www.ncedc.org/). 


40°  -4 


Figure  32:  Map  of  the  western  United  States.  Red  stars  are  Nevada  Test  Site  (NTS) 
explosions,  black  circles  are  earthquakes  and  the  white  triangles  are  broadband  seismic 
stations. 
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3.4.3  Inversion  Procedure 


The  MT  inversion  follows  the  method  developed  by  Minson  and  Dreger  (2008) 
where  the  complete  waveform  is  inverted  in  a  weighted  linear  least  square  formulation.  The 
displacement  waveform  is  expressed  as  a  linear  combination  of  the  GFs  weighted  by  the  MT 
elements.  The  first  step  is  data  preparation  where  the  hypocenter  and  station  location  are 
specified,  and  the  data  are  processed  as  described  in  the  beginning  of  this  section.  After 
reading  the  data  files,  the  MT  inversion  algorithm  then  loads  the  pre-computed  GFs  for  each 
station,  and  performs  a  waveform  cross-correlation  between  the  data  and  GF  at  each  station 
to  determine  the  initial  time  shifts  (Pasyanos  et  ah,  1996).  A  more  detailed  description  on 
computing  the  GF  is  provided  in  the  next  subsection.  The  time  shift  that  gives  the  highest 
correlation  between  data  and  GFs  is  selected  (Pasyanos  et  al.,  1996).  We  do  not  allow  for 
time  shifts  greater  than  half  a  cycle  of  the  minimum  wavelength  to  avoid  cycle  shifting. 
Applying  time  shifts  prior  to  the  actual  inversion  is  a  common  procedure  used  to  account  for 
imperfect  travel  time  predicted  by  the  model  and  errors  in  event  origin  time.  The  final  step  is 
the  inversion  where  we  do  not  constrain  the  inversion  to  have  zero  trace,  instead  it  solves  for 
all  six  components  of  the  seismic  MT. 


3.4.4  Reciprocal  Green’s  Functions 

We  compute  elastodynamic  GFs  for  each  source-receiver  pair  using  ID  and  3D  Earth 
models.  The  GFs  are  the  impulse  response  whose  orientation  is  given  by  the  basic  MT 
elements  or  by  the  fundamental  fault  types  (e.g.  Helmberger,  1983;  Jost  and  Herrmann, 
1989).  The  three  basic  faults  types  are  the  vertical  dip-slip  (DS),  45-degree  dip-slip  (DD), 
and  vertical  strike-slip  (SS).  Any  arbitrary  (full)  MT  can  be  described  as  a  linear 
combination  of  the  six  basic  MT  elements  ( Mxx ,  Mxy,  Mxz,  Myy,  Myz,  Mzz),  where  the 
first  index  represents  the  direction  of  the  force,  and  the  second  the  direction  of  a  spatial 
derivative.  Thus  Mxx,  Myy  and  Mzz  represent  the  strength  of  vector  dipole  forces,  and 
Mxy,  Mxz,  Myz,  Myx,  Mzx  and  Mzy  represent  the  strength  of  force  couples. 
Conservation  of  angular  momentum  leads  to  symmetry  of  the  MTs,  where  Mxy  =  Myx , 
yielding  the  6  independent  elements.  For  the  deviatoric  MT,  where  the  trace  is  zero,  five 
basic  tensor  elements  or  the  three  fundamental  fault  types  can  describe  the  seismic  tensor 
(e.g.  Langston,  1981;  Langston  and  Helmberger,  1975). 


The  GFs  are  computed  for  each  source-receiver  pair  and  processed  the  same  way  as 
the  data.  The  ID  GFs  are  computed  using  the  western  US  velocity  model  (Song  et  al., 
1996)  and  frequency  wavenumber  integration  method  (Wang  and  Herrmann,  1980; 
Herrmann  and  Wang,  1985;  Herrmann,  2013).  The  3D  GFs  are  computed  using  the 
anelastic  FD  code,  SW4,  developed  at  the  Center  of  Applied  Scientific  Computing  at 
Lawrence  Livermore  National  Laboratory  (LLNL;  part  of  the  Computational  Infrastructure 
for  Geodynamics).  The  conventional  approach  to  compute  3D  GFs  using  FD  is 
computationally  expensive  because  we  need  to  simulate  wave  propagation  at  regional 
distances  with  attenuation  and  at  shorter  periods.  For  each  receiver,  the  number  of 
simulations  required  to  compute  3D  GFs  is  six  times  the  number  of  sources  to  compute  the 
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six  basic  MT  elements.  Therefore  to  reduce  the  computation  cost  we  apply  the  reciprocity 
principle  to  compute  the  3D  GFs.  Source-receiver  reciprocity  for  seismic  waves  is  given  by 
the  Betti’s  theorem  (Aki  and  Richards,  2002;  Dahlen  and  Tromp,  1998)  which  states  that 
the  location  of  the  source  and  observation  can  be  switched  but  still  produce  the  same  exact 
elastic  response.  Eisner  and  Clayton  (2001)  and  Graves  and  Wald  (2001)  implemented 
source-receiver  reciprocity  using  the  FD  approach  to  simulate  ground  motion  for  any 
number  of  arbitrary  sources  from  a  set  of  stations.  The  advantage  of  using  reciprocity  rather 
than  the  forward  approach  is  that  it  is  only  necessary  to  run  three  calculations  at  every 
receiver  location.  Commonly  there  are  fewer  stations  than  sources  in  a  typical  monitoring 
scenario.  We  compute  the  elastic  response  at  the  source  due  to  a  single  force  in  three 
orthogonal  directions  (X,  Y  and  Z)  at  the  receiver.  For  each  reciprocal  single  force  the 
numerical  spatial  derivative  of  the  displacement  in  the  interior  of  the  model  can  be 
computed.  Then  it  is  possible  to  combine  the  derivatives  to  construct  the  individual  MT 
components.  For  a  single  force  oriented  along  the  X  direction: 


G 


x 

XX 


dy 


(27) 


du y  dux 

— — - 

dx  dy 


dux  dux 
— —  +  — - 

dx  dz 


dy 


dif 

_ y_ 

dz 


There  are  similar  equations  for  unit  forces  oriented  in  the  Y  and  Z  directions.  We 
call  the  elastic  response  computed  via  reciprocity  the  reciprocal  Green’s  functions 
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(RGF)  to  distinguish  them  from  the  conventional  (forward)  approach.  We  can  then 
compute  the  X  component  of  displacement  for  an  arbitrary  source  as  the  weighted  sum  of 
the  RGFs: 


ux  =M  Gj  +M  Gx  +M  Gx  +M  Gx  +M  Gx  +M  Gx 

xx  xx  yy  yy  zz  zz  xy  xy  xz  xz  yz  yz 


(28) 


The  capitol  superscripts  refer  to  the  component  of  motion,  which  corresponds  to  the 
direction  of  the  applied  unit  force  at  the  receiver  position.  Given  the  strike,  rake,  dip  and 
scalar  moment  of  the  fundamental  fault  types  we  can  compute  the  /\z/,y  coefficients  (Eq.  1  in 
Box  4.4  of  Aki  and  Richards,  2002),  and  substituting  them  into  Eq.  5.2  we  can  construct 
the  displacement  response  for  each  fundamental  fault  type  Green’s  function  for  the 
moment  tensor  inversion. 

The  3D  isotropic  model  by  Moschetti  et  al.  (2010)  is  linearly  interpolated  in  the 
vertical  direction  and  smoothing  is  applied  in  the  horizontal  direction  using  Gaussian 
averaging.  The  smoothing  is  necessary  because  the  velocity  model  is  much  coarser  than 
the  computation  grid.  For  the  reciprocal  FD  calculations  we  use  a  grid  spacing  of  250  m 
and  a  Gaussian  source  function  with  a  duration  of  0.5  seconds.  We  did  not  include  any 
water  layers  or  topography  in  the  3D  calculations,  and  implemented  attenuation  using  a 
function  of  local  Vs  [lcm/s]  (e.g.  Graves  et  al.,  2008;  Olsen  et  al.,  2009),  such  that  Qs  = 
50  Vs  and  Qp  =  2Qs.  Figure  33  shows  a  comparison  between  a  forward  FD  calculation  and 
a  reciprocal  calculation.  In  this  case,  the  source  depth  is  1  km  and  the  receiver  is  100 
km  away  and  at  the  surface.  There  is  excellent  agreement  between  the  forward  and 
reciprocal  simulations. 


The  FD  calculations  were  performed  on  the  Fivermore  Computing  Center  (FC)  8- 
core  Xeon  35-2670  Finux  cluster,  cab,  consisting  of  1,296  nodes,  each  with  16  cores  per 
node  and  32  GB  memory.  Ideally  we  should  use  the  same  computation  domain  for  all  the 
reciprocal  calculations  but  the  run  time  increases  rather  quickly  at  large  distances  due  to  a 
drastic  increase  in  the  number  of  grid  points  and  the  duration  of  the  simulation.  Because 
the  3D  model  is  smooth  and  we  are  only  comparing  waveforms  as  short  as  8  seconds,  we 
justify  using  different  computation  domains  by  positioning  the  source  and  receivers  to  be  at 
least  250  grid  points  away  from  the  boundary  and  use  a  sufficiently  large  absorbing 
boundary  layer  (30  km  in  all  directions)  to  minimize  boundary  reflections.  Running  on  512 
processors,  the  run  time  for  one  RGF  calculation  with  approximately  1 .2  billion  grid  points 
and  1,1924  time  steps  takes  about  10  hours. 
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Figure  33:  A  comparison  between  forward  and  reciprocal  finite- difference  calculations  at 
100  km  distance.  The  wavefonns  are  in  velocity  (cm/s)  and  bandpass-filtered  between  0.02 
and  0.15  Hz  with  a  4-pole  acausal  Butterworth  filter.  The  black  lines  are  the  reciprocal 
calculations  and  the  brown  dashed  lines  are  the  conventional  (forward)  calculations. 


3.4.5  Moment  Tensor  Inversion  Results 

We  perform  ID  and  3D  MT  inversion  and  compute  Network  Sensitivity  Solutions 
[NSS]  (e.g.  Ford  et  al.,  2010)  for  all  events.  For  four  out  of  the  six  events,  we  also  compare 
the  source  properties  in  two  different  frequency  bands.  For  the  earthquakes  (Armagosa  and 
Little  Skull  Mountain),  we  search  for  the  best  depth  using  the  ID  model  and  invert  for  ID 
and  3D  MT  solutions  at  these  depths,  whereas  for  the  explosions  (METROPOLIS,  COSO, 
HOYA  and  JUNCTION),  we  fix  the  source  depth  at  1  km  for  both  ID  and  3D  GFs.  For  the 
comparison  at  different  frequency  bands,  we  use  data  from  the  BDSN  and  SCSN  networks, 
and  when  available,  we  also  compute  solutions  including  the  LLNL  network  stations  but 
since  the  raw  data  and  instrument  responses  for  these  stations  are  currently  not  open  we 
cannot  compare  the  solutions  in  the  different  frequency  bands. 

Low  Signal-to-Noise  Events 

COSO  and  Amargosa  have  relatively  low  SNR  compare  to  the  other  events  in  this 
study,  therefore  we  did  not  compare  the  ID  and  3D  inversion  results  at  different  frequency 
bands,  instead  the  inversion  is  done  between  10  to  50  seconds.  Figure  34a  shows  the  full 
MT  inversion  results  for  COSO.  The  ID  solution  fits  the  data  slightly  better  than  the  3D 
solution,  while  the  moment  magnitude  ( Mw )  using  the  3D  model  is  slightly  higher  than  that 
of  the  ID  model.  Both  the  ID  and  3D  GF  inversions  recover  solutions  that  are  dominated 
by  the  explosion  component,  representing  65%  and  68%  of  the  total  seismic  moment, 
respectively.  On  the  other  hand,  the  off-diagonal  components  of  the  MTs  are  quite  different, 
where  the  ID  solution  has  a  deviatoric  component  that  is  mostly  a  normal  DC  mechanism 
and  the  3D  solution  is  a  mixture  of  vertical  DS  and  CLVD  mechanisms.  Thus  while  the 
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explosive  nature  of  the  event  is  clearly  recovered  in  both  cases  there  is  velocity  model 
dependence  on  the  lesser  “tectonic”  release  component  of  the  source. 

For  the  Amargosa  earthquake  (Figure  34b)  the  ID  solution  is  fitting  17%  better  than 
the  3D  solution,  a  much  greater  difference  compared  to  COSO.  This  is  interesting  in  of 
itself  because  it  clearly  shows  that  the  use  of  a  3D  velocity  model,  including  one  that  is 
constrained  by  data  (surface  wave  dispersion)  will  not  necessarily  lead  to  improved  ability  to 
fit  the  data.  The  ID  and  3D  fault  plane  orientations  are  also  quite  different  where  the  ID 
solution  is  closer  to  a  normal  mechanism.  In  contrast,  the  3D  solution  has  a  vertical  DS 
component  that  is  significantly  rotated  with  respect  to  the  ID  solution.  In  addition,  the  3D 
solution  results  in  a  relatively  large  “false”  ISO  of  36%  of  the  total  seismic  moment. 
However,  the  ISO  component  is  not  statistically  significant  as  determined  by  the  F-test 
(Menke,  1989;  Menke,  2012;  Templeton  and  Dreger,  2006),  where  the  level  of  statistical 
significance  with  the  additional  degree  of  freedom  in  the  full  MT  inversion  is  only  50%. 
Nevertheless,  this  illustrates  that  imperfect  3D  velocity  models  can  also  increase  the  ISO 
moment  in  earthquakes,  contrary  to  the  results  of  Covellone  and  Savage  (2012).  Although 
we  see  an  increase  in  ISO  component  for  short  period  inversions  the  ID  (Figure  35a)  and  3D 
NSS  (Figure  35b)  exhibit  a  typical  earthquake-like  signature  with  a  bullseye  pattern  (Ford  et 
al.,  2010). 

Little  Skull  Mountain  Earthquake 

The  Md  4.3  2002  Little  Skull  Mountain  earthquake  occurred  near  NTS, 
approximately  6  km  southwest  of  the  1992  Little  Skull  Mountain  mainshock.  At  20  to  50 

seconds,  we  have  data  from  five  stations  with  an  azimuthal  coverage  of  89  .  For  the 
comparison  at  long  period  (Figure  36a-c),  both  ID  and  3D  MT  solutions  fit  the  data  equally 
well,  and  have  similar  mechanisms.  The  solution  is  a  normal  mechanism  and  is 
predominantly  DC.  The  DC  component  for  the  3D  solution  is  slightly  higher  than  the  ID 
solution,  and  unlike  the  explosions  the  ID  solution  Mw  is  higher  than  the  3D  solution. 

For  the  inversion  at  short  period  we  did  not  include  MHC  due  to  low  SNR.  Using 
data  from  8  to  20  second  yields  ID  and  3D  MTs  with  only  50%  and  16%  DC,  respectively 
(Figure  36a)  considerably  lower  than  solutions  from  the  20  to  50  second  passband.  The  low 
DC  in  the  3D  solution  is  particularly  alarming.  The  3D  solution  also  recovers  a  rotated 
normal  mechanism  with  respect  to  long  period  MTs  (about  40  degrees).  The  Mw  is  similar 
for  all  solutions  and  the  ID  model  performs  better  than  the  3D  model  mainly  because  the 
3D  model  cannot  fit  the  tangential  component  as  well  as  the  ID  model  (Figure  36d-e).  We 
also  see  an  increase  in  variance  at  CMB  using  the  3D  model.  We  identify  ISA  to  be  a 
problematic  station  in  the  3D  inversion  because  the  synthetics  cannot  fit  the  data  well  with 
time  shifts  less  than  half  the  minimum  wavelength  (4s).  Therefore  the  solution  with  the 
highest  VR  within  the  4-second  window  is  a  solution  with  zero  shifting.  Considering  a 
three-station  inversion  without  ISA  (Figure  37)  we  obtain  a  solution  with  52%  DC  and  a 
normal  mechanism  consistent  with  long  period  MTs  and  the  ID  short  period  MT. 
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Figure  34:  COSO  and  Amargosa  full  moment  tensor  inversion  results  with  ID  and  3D 
Green’s  functions  filtered  between  10  to  50  seconds.  Full  moment  tensor  focal  mechanism 
and  the  deviatoric  component  of  the  solution  are  plotted  as  well  as  the  moment  magni¬ 
tude  (Mw),  percent  double-couple  (DC),  percent  compensated  linear  vector  dipole 
(CLVD),  percent  isotropic  (ISO)  and  variance  reduction  (VR).  VR  from  deviatoric 
inversion  is  in  parentheses. 


Whether  or  not  we  include  ISA  in  the  3D  MT  inversion,  the  ISO  component  remains 
high,  up  to  40%,  which  is  quite  large  for  an  earthquake.  However  the  F-test  reveals  the 
ISO  component  is  not  statistically  significant  with  a  50%  confidence  level  of  significance 
for  inversions  with  and  without  ISA.  In  addition,  ID  (Figure  38a)  and  3D  NSS  (Figure 
38b-c)  all  exhibit  an  earthquake-like  distribution  and  similar  level  of  solution  uncertainty 
where  the  normalized  VR>95%  (67%,  43%  and  54%,  respectively)  are  about  the  same 
size.  Similar  to  Amargosa,  the  results  show  that  imperfect  3D  velocity  models  can  also 
increase  the  ISO  moment  in  earthquakes. 


Larger  NNSS  (NTS)  Explosions 

We  compare  ID  and  3D  MT  inversions  for  three  NTS  explosions  METROPOLIS, 
HOYA  and  JUNCTION.  Of  the  three  explosions  HOYA  has  strong  SH  radiation  on  the 
tangential  component,  indicative  of  a  larger  tectonic  release.  For  METROPOLIS  we  have 
data  from  three  stations  PAS,  CMB  and  MHC  where  we  can  compare  the  solutions  at 
different  frequency  bands.  Additionally,  we  are  able  to  compute  a  solution  including  two 
additional  stations,  MNV  and  LAC  from  the  LLNL  network,  increasing  the  azimuthal 
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coverage  from  75  to  123  .  Data  from  MNV  and  LAC  are  filtered  between  10  to  30 
seconds.  PAS  and  CMB  have  high  SNR  at  both  long  and  short  periods  whereas  MHC  has 
relatively  noisy  data  at  long  periods  for  the  two  horizontal  components,  but  relatively  high 
SNR  for  the  vertical  component.  At  20  to  50  seconds,  we  see  very  little  SH  energy  in  the 
tangential  component  and  essentially  none  at  CMB.  The  deviatoric  part  of  the  solution  is 
very  similar  for  both  ID  and  3D  MTs;  however,  the  3D  solution  has  a  much  higher  ISO 
component,  16%  higher  than  the  ID  solution  (Figure  39a).  Again,  the  3D  solution  has  a 
higher  Mw  that  is  closer  to  the  reported  magnitude  (Md)  in  the  catalog  (Table  5.1).  Although 
the  VR  improved  slightly  for  all  three  stations  from  ID  to  3D  (Figure  39b-c)  the  increase  is 
only  a  few  percent.  At  8  to  20  seconds,  the  ISO  component  of  the  3D  solution  decreased  to 
59%,  in  contrast  the  ISO  component  of  the  ID  solution,  which  increased  to  81%.  The 
deviatoric  part  of  the  solution  changed  significantly  from  long  period  to  short  period,  instead 
of  a  reverse  sense  of  motion  for  both  solutions  at  long  periods,  at  short  periods  the  ID 
solution  is  mostly  DS  and  the  3D  solution  has  a  large  CLVD  component.  In  terms  of 
waveform  fits  (Figure  39d-3)  overall  the  ID  model  is  fitting  the  data  better  than  the  3D 
model  due  to  significant  increase  in  VR  at  CMB.  On  the  other  hand,  the  3D  model  performs 
better  along  the  path  to  MHC.  Mw  decreased  for  both  solutions  from  the  long  period  case  to 
short  period  case,  but  again  the  3D  solution  has  a  higher  Mw.  The  best  solution  is  the 
inversion  using  data  from  five  stations  and  filtered  between  10  to  50  seconds  (Figure  39a). 
Both  ID  and  3D  results  have  high  ISO  components  and  similar  mechanisms  for  the  off- 
diagonal  components.  Mw  for  the  3D  solution  is  0.3  magnitude  units  higher  than  the  ID 
solution,  and  overall  the  3D  model  is  fitting  a  few  percent  better  than  the  ID  model. 


a.  1DNSS,  10  -50  s 


b.  3D  NSS,  10 -50  s 


Figure  35:  Amargosa  NSS  for  (a)  ID  and  (b)  3D  Green’s  functions  filtered  between  10 
to  50  seconds.  White  circle  marks  the  location  of  best  full  moment  tensor  in  source- 
type  space. 
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a.  2002  Little  Skull  Mountain  Moment  Tensor  Solutions 
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Figure  36:  2002  Little  Skull  earthquake  full  moment  tensor  inversion  results  with  ID  and  3D 
Green’s  functions  and  at  different  frequency  bands,  (a)  Full  moment  tensor  focal  mechanism 
and  the  deviatoric  component  of  the  solution  are  plotted  as  well  as  the  moment  magnitude 
(Mw),  percent  double-couple  (DC),  percent  compensated  linear  vector  dipole  (CLVD), 
percent  isotropic  (ISO)  and  variance  reduction  (VR).  VR  from  deviatoric  inversion  is  in 
parentheses,  (b-e)  Data  (solid  line)  and  synthetic  waveforms  (dashed  line)  plotted  from  left 
to  the  right  are  the  tangential,  radial  and  vertical  components. 
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a.  2002  Little  Skull  Mountain  Moment  Tensor  Solutions  (3-station) 
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Figure  37:  2002  Little  Skull  earthquake  three-station  inversion  using  8  to  20  second  3D 
Green’s  functions,  (a)  Full  moment  tensor  focal  mechanism  and  the  deviatoric  component 
of  the  solution  are  plotted  as  well  as  the  moment  magnitude  ( Mw ),  percent  double-couple 
(DC),  percent  compensated  linear  vector  dipole  (CLVD),  percent  isotropic  (ISO)  and 
variance  reduction  (VR).  VR  from  deviatoric  inversion  is  in  parentheses,  (b)  Data  (solid 
line)  and  synthetic  waveforms  (dashed  line)  plotted  from  left  to  the  right  are  the 
tangential,  radial  and  vertical  components. 


Unlike  the  other  explosions  in  this  study,  HOYA  has  strong  Love  wave  energy 
across  all  stations.  The  data  also  have  high  SNR  in  the  frequency  band  we  examined 
and  better  station  coverage  compare  to  the  other  explosions.  The  azimuthal  coverage  for 
HOYA  is  100  °  for  the  five-station  inversion,  and  the  best  solution  including  two  additional 
LLNL  stations  MNV  and  LAC  has  an  azimuthal  coverage  of  132  ° .  For  solutions  from  20  to 
50  seconds  (Figure  40a-c),  the  ISO  component  is  around  60%  for  both  velocity  models; 
however  the  3D  solution  has  a  larger  DC  component  and  a  higher  Mw.  The  non-ISO 
mechanism  is  quite  different  between  the  two  models,  the  ID  case  is  a  normal  mechanism 
and  the  3D  case  is  a  vertical  DS  mechanism.  The  3D  model  fits  the  data  better  than  the  ID 
model.  As  seen  in  Figure  40b-c  the  amplitudes  between  the  data  and  synthetics  at  stations 
PAS,  PFO  and  BKS  are  in  better  agreement  when  we  use  a  3D  model. 
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a.  1 D  NSS,  8  -  20  s,  4-station  b.  3D  NSS,  8  -  20  s,  4-station 


c.  3D  NSS,  8  -  20  s,  3-station 


Figure  38:  2002  Little  Skull  Mountain  NSS  using  8  to  20  second  Green’s  functions  for  (a)  ID 
4-station  inversion,  (b)  3D  4-station  inversion  and  (c)  3D  3-station  inversion.  White  circles 
mark  the  location  of  best  full  moment  tensor  in  source-type  space. 


From  8  to  20  seconds  (Figure  40a, d-e)  the  3D  solution  has  an  ISO  component  of  73% 
whereas  the  ID  solution  has  an  ISO  component  of  54%.  The  3D  model  cannot  fit  the  data 
well  on  the  tangential  components  at  all  stations  except  for  BKS.  We  can  see  the  waveform 
fits  are  significantly  better  along  the  path  to  BKS  when  using  3D  GFs  instead  of  ID  GFs, 
resulting  in  higher  overall  VR  for  the  3D  solution.  At  short  periods  the  deviatoric  part  of 
the  solution  has  a  different  mechanism  compared  to  the  long  period  solutions.  Although  the 
3D  case  still  consists  of  a  vertical  DS  mechanism  the  fault  orientation  has  rotated  about 

45  °  towards  the  southwest.  In  contrast  to  the  other  explosions  the  Mw  is  slightly  higher 
when  using  ID  GFs  at  short  periods.  The  long  period  solutions  are  more  consistent  with 
the  seven-station  inversion  because  the  frequency  band  used  for  the  best  solution  is  between 
10  to  50  seconds.  We  use  inverse  variance  weighting  instead  of  inverse  distance  weighting 
for  the  seven-station  inversion  because  data  from  MNV  and  LAC  are  filtered  between  10  to 
30  seconds,  resulting  in  higher  amplitudes,  while  the  rest  of  the  stations  are  filtered  between 
20  to  50  seconds.  For  the  seven-station  inversion,  the  3D  solution  has  a  higher  Mw  and  ISO 
component,  and  a  vertical  DS  mechanism  instead  of  a  normal  mechanism  for  the  deviatoric 
part  of  the  full  MT. 


For  JUNCTION,  we  have  data  from  four  stations  BKS,  PAS,  PFO  and  ISA  for  the 
comparison  at  different  periods,  and  there  is  an  additional  LLNL  station  LAC  in  the  five- 
station  inversion.  For  the  five-station  inversion,  all  data  from  the  BDSN  and  SCSN  network 
stations  are  filtered  between  10  to  50  seconds  and  LAC  is  filtered  between  10  to  30  seconds. 
Comparing  the  results  at  20  to  50  seconds  (Figure  41a)  full  MT  solutions  using  either  model 
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have  excellent  waveform  fits;  the  3D  solution  has  a  VR  of  95%  and  the  ID  model  has  a  VR 
of  83%.  VR  increases  for  all  stations  using  3D  GFs  but  the  biggest  improvements  are  along 
the  paths  to  BKS  and  PFO  (Figure  41b-c).  Similar  to  HOYA,  the  amplitudes  between  data 
and  synthetics  on  the  radial  and  vertical  components  agree  well  when  3D  GFs  are  used  but 
not  as  well  on  the  tangential  component  for  BKS  and  PFO.  Both  solutions  are 
predominantly  ISO  but  the  3D  solution  has  an  ISO  component  17%  higher  than  the  ID 
solution.  Mw  from  the  3D  solution  is  also  higher  than  the  ID  case.  The  deviatoric 
component  of  the  solutions  have  different  fault  orientations  and  sense  of  motion  between  the 
ID  and  3D  case;  the  ID  solution  consists  of  a  normal  mechanism  and  the  3D  solution 
consists  of  aCLVD  mechanism  with  the  major  vector  dipole  in  tension  (+CLVD). 

Comparing  the  solutions  at  periods  between  8  and  20  seconds,  we  see  again  that  the 
paths  to  BKS  and  PFO  are  well  modeled  using  3D  GFs,  which  are  consistent  with  the  results 
at  long  periods.  However,  the  VR  decreases  for  the  paths  to  ISA  and  PAS  when  using  3D 
GFs  (Figure  41d-e).  Different  than  the  results  at  long  period,  the  ID  solution  has  a  larger 
ISO  component  of  75%  while  the  3D  solution  has  an  ISO  component  of  69%.  The  non-ISO 
part  of  the  MT  consists  of  a  DS  mechanism  but  with  different  orientations  between  ID  and 
3D  cases.  Jackknife  resampling  of  different  station  combinations  shows  PAS  is  a  key  station 
in  the  inversion  at  short  period,  the  ISO  component  decreases  when  we  remove  PAS, 
especially  for  the  ID  case  where  the  ISO  disappears  and  the  solution  is  a  CLVD 

mechanism  with  the  major  vector  dipole  in  compression  (—CLVD).  The  best  solution  using 

data  from  five  stations  are  mostly  ISO  but  for  the  ID  case  the  ISO  component  is  reduced  to 
56%.  In  terms  of  Mw,  the  3D  solution  overestimates  the  magnitude  when  we  compare  it 
to  the  reported  ML  in  Springer  et  al.  (2002).  The  deviatoric  part  of  the  solution  comprises 
of  a  vertical  DS  mechanism  for  the  ID  case  and  a  +CLVD  mechanism  for  the  3D  case.  The 
drop  in  VR  for  the  ID  solution  is  due  to  lower  fits  at  station  LAC  whereas  all  other  stations 
have  VR  above  70%. 
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a.  METROPOLIS  Moment  Tensor  Solutions 
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Figure  39:  METROPOLIS  full  moment  tensor  inversion  results  with  ID  and  3D  Green’s 
functions  and  at  different  frequency  bands,  (a)  Full  moment  tensor  focal  mechanism  and  the 
deviatoric  component  of  the  solution  are  plotted  as  well  as  the  moment  magnitude  ( Mw ), 
percent  double-couple  (DC),  percent  compensated  linear  vector  dipole  (CLVD),  percent 
isotropic  (ISO)  and  variance  reduction  (VR).  VR  from  deviatoric  inversion  is  in  parentheses, 
(b-e)  Data  (solid  line)  and  synthetic  waveforms  (dashed  line)  plotted  from  left  to  the  right 
are  the  tangential,  radial  and  vertical  components. 
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a.  HOYA  Moment  Tensor  Solutions 
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Figure  40:  HOYA  full  moment  tensor  inversion  results  with  ID  and  3D  Green’s  functions 
and  at  different  frequency  bands,  (a)  Full  moment  tensor  focal  mechanism  and  the 
deviatoric  component  of  the  solution  are  plotted  as  well  as  the  moment  magnitude  (Mw), 
percent  double-couple  (DC),  percent  compensated  linear  vector  dipole  (CLVD),  percent 
isotropic  (ISO)  and  variance  reduction  (VR).  VR  from  deviatoric  inversion  is  in 
parentheses,  (b-e)  Data  (solid  line)  and  synthetic  wavefonns  (dashed  line)  plotted  from  left 
to  the  right  are  the  tangential,  radial  and  vertical  components. 
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a.  JUNCTION  Moment  Tensor  Solutions 
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Figure  41:  JUNCTION  full  moment  tensor  inversion  results  with  ID  and  3D  Green’s 
functions  and  at  different  frequency  bands,  (a)  Full  moment  tensor  focal  mechanism  and 
the  deviatoric  component  of  the  solution  are  plotted  as  well  as  the  moment  magnitude 
( Mw ),  percent  double-couple  (DC),  percent  compensated  linear  vector  dipole  (CLVD), 
percent  isotropic  (ISO)  and  variance  reduction  (VR).  VR  from  deviatoric  inversion  is  in 
parentheses,  (b-e)  Data  (solid  line)  and  synthetic  wavefonns  (dashed  line)  plotted  from  left 
to  the  right  are  the  tangential,  radial  and  vertical  components. 
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Comparison  in  Source-Type  Space 

Network  Sensitivity  Solutions  (NSS)  utilizing  source-type  representations  proposed 
by  Hudson  et  al.  (1989)  and  Tape  and  Tape  (2012a)  and  Tape  and  Tape  (2012b)  have  been 
developed  and  used  in  previous  studies  to  estimate  uncertainties  in  source  inversions  (e.g. 
Ford  et  al.,  2010;  Guilhem  et  al.,  2014;  Chiang  et  al.,  2014;  Nayak  and  Dreger,  2015). 
We  would  like  to  compare  the  estimated  uncertainties  in  the  source-type  space  between  ID 
and  3D  GFs  to  see  if  the  3D  model  reduces  the  source-type  uncertainties,  therefore 
providing  better  constraints  on  source  mechanisms  and  increasing  confidence  in 
discrimination.  We  compute  the  NSS  for  Little  Skull  Mountain,  METROPOLIS,  HOYA 
and  JUNCTION  using  the  same  station  configuration  and  filter  parameters  in  the  MT 
inversion  analysis.  The  results  comparing  ID  and  3D  NSS  at  two  frequency  bands  are 

presented  in  Figure  42,  where  we  plotted  solutions  with  normalized  VR>95%.  The 

comparison  for  the  2002  Little  Skull  Mountain  earthquake  shows  at  long  periods  (20-50 
seconds)  the  spread  of  the  MT  uncertainties  are  larger  using  3D  GFs,  but  at  short  periods 
(8-20  seconds)  the  spread  is  about  the  same  between  the  two  models.  If  we  look  at  the 
distribution  of  the  source-type  parameters  y  and  5  (Figure  43a)  we  see  a  wider 
distribution  in  the  volumetric  component  (5)  of  the  MT  when  inverting  using  3D  GFs.  At 
higher  frequencies,  the  two  NSS  populate  different  areas  of  the  source-type  space;  the  3D 
NSS  extends  into  the  region  of  positive  volume  change  while  the  ID  NSS  extends  into  the 
region  of  negative  volume  change.  There  are  noticeable  shifts  in  the  mean  of  the 
distributions  for  both  ID  and  3D  NSS  (Figure  43a),  where  the  distributions  at  short 
periods  deviate  away  from  a  pure  DC  mechanism.  For  the  explosions,  ID  and  3D  NSS 
(Figure  42)  at  long  periods  show  the  spread  of  the  estimated  uncertainties  associated  with 
source  inversion  varies  between  events:  the  uncertainties  for  JUNCTION  are  similar  using 
either  ID  or  3D  GFs,  for  HOYA  the  solutions  are  better  constrained  using  3D  GFs  but  the 
results  are  opposite  for  METROPOLIS.  At  short  periods  we  see  little  difference  between 
the  ID  and  3D  NSS,  except  for  HOYA  where  the  solution  is  very  well-constrained  using 
ID  GFs.  The  distribution  of  y  and  5  for  explosions  are  negative  skewed  (Figure  43b) 
meaning  they  have  a  long  tail  in  the  negative  direction,  for  both  ID  and  3D  models  and  at 
both  long  and  short  periods,  the  exception  is  the  HOYA  ID  NSS  where  y  is  positive 
skewed  resulting  in  a  population  of  solutions  near  the  theoretical  explosion/opening  crack 

mechanisms.  All  explosion  NSS  have  a  mean  5  around  30  °  and  we  do  not  see  a  large 
difference  in  the  mean  5  values  between  ID  and  3D  NSS. 
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Figure  42.  Comparison  of  Network  Sensitivity  Solution  (NSS)  of  four  different  events:  2002 
Little  Skull  earthquake,  METROPOLIS,  HOYA  and  JUNCTION.  The  shaded  regions  and 
contour  lines  show  the  populations  of  solutions  with  normalized  VR>95%  for  long  period 
and  short  period  waveform  inversions,  respectively.  The  blue  and  pink  colors  are  solutions 
computed  using  3D  and  ID  Green’s  functions,  respectively. 
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a.  2002  Little  Skull  Earthquake 
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Figure  43:  Comparisons  of  source-type  parameters  y  and  5  for  two  events  (a) 
2002  Little  Skull  Mountain  earthquake  and  (b)  NTS  explosion  JUNCTION. 
Histograms  show  y  and  5  from  waveform  inversions  at  two  frequency  bands, 
the  gray  bars  represent  the  use  of  ID  Green’s  functions,  the  white  bars 
represent  the  use  of  3D  Green’s  functions,  and  the  gray  and  black  lines  are  the 
mean  values  from  ID  and  3D  inversions,  respectively. 
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3.4.6  Discussion 


Comparisons  between  ID  and  3D  source  inversions  show  differences  in  source 
mechanism  in  both  the  long  period  and  short  period  cases.  For  the  20  to  50  second  case, 
the  ID  and  3D  solutions  for  the  2002  Little  Skull  Mountain  earthquake  show  minimal 
difference  in  source  mechanism  and  magnitude.  In  contrast,  the  explosions  show  an 
increase  in  magnitude  and  waveform  fits  from  ID  to  3D  inversions.  The  off-diagonal 
components  are  not  well-constrained  and  vary  depending  on  the  frequency  band  and 
velocity  model  used.  Contrary  to  Hingee  et  al.  (2011)  and  Covellone  and  Savage  (2012), 
where  MT  solutions  were  compared  from  40  to  200  seconds  and  25  to  125  seconds, 
respectively,  in  general  we  see  a  reduction  in  waveform  fits  for  the  3D  inversions  at 
relatively  short  periods  (8  to  20  seconds)  and  for  the  two  earthquakes  we  see  an  increase  in 
non-DC  components,  suggesting  finer  details  of  the  Earth’s  structure  may  not  be  well- 
represented  by  the  3D  model.  The  exceptions  are  the  paths  to  BKS  and  MHC,  particularly 
BKS,  where  we  see  significant  improvements  in  waveform  fits  when  3D  GFs  are  used  in 
the  short  period  inversions.  Paths  crossing  the  extensional  regimens  and  low  velocity 
sediments  (Figure  44a)  are  better  represented  by  the  3D  model,  where  stations  BKS  and 
MHC  are  situated  on  top  of  the  low  velocity  zone  (Figure  44b).  Based  on  the  four  events 
analyzed  in  this  study,  the  long  period  time  domain  MT  inversion  results  in  higher 
wavefonn  fits  when  3D  GFs  are  applied,  however  we  do  not  see  a  significant  reduction  in 
uncertainties  associated  with  the  MT  using  synthetics  derived  from  the  3D  model. 
JUNCTION  has  little  difference  in  source  uncertainties  between  the  ID  and  3D  case,  Little 
Skull  Mountain  and  METROPOLIS  show  higher  uncertainties  for  the  3D  case  whereas 
HOYA  shows  lower  uncertainties  for  the  3D  case.  A  larger  sample  size  is  required  to  make 
more  useful  interpretations  about  the  use  of  3D  models  in  estimating  MT  source 
uncertainties. 
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Figure  44:  (a)  Vs  at  10  km  depth  where  surface  waves  are  most  sensitive  to  the  structure 
at  these  period  ranges,  (b)  A  cross-sectional  view  of  crustal  and  upper  mantle  velocities 
across  A- A’. 

Moschetti  et  al.  (2010)  noted  the  surface  wave  model  has  larger  uncertainties  in 
shear  wave  speeds  in  extensional  regions  across  the  western  US,  near  the  Moho  (lower 
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crust)  and  the  shallowest  parts  of  the  crust.  There  is  a  significant  trade-off  between  crustal 
thickness  and  shear  wave  speed  resulting  in  an  increase  in  shear  wave  uncertainties  around 
the  Moho  (between  35  to  45  km  depth).  In  our  3D  short  period  inversions,  Rayleigh  waves 
are  fitting  better  than  the  Love  waves.  The  3D  model  used  was  determined  using  dispersion 
data  from  Rayleigh  wave  group  and  phase  velocities  and  Love  wave  phase  velocities,  Love 
wave  group  velocities  are  not  included  due  to  large  uncertainties.  Also,  the  Love  wave 
dispersion  maps  are  in  period  bands  of  8  to  32  seconds  whereas  the  Rayleigh  wave 
dispersion  maps  are  in  period  bands  of  6  to  40  seconds.  Therefore  in  addition  to  a  priori 
constraints  on  sediment  and  crustal  thicknesses  (Laske  and  Masters,  1997;  Gilbert  and 
Fouch,  2007)  the  details  of  the  3D  model  were  constrained  mostly  by  short  period 
Rayleigh  wave  data.  In  comparison,  Hingee  et  al.  (2011)  used  a  3D  radially  anisotropic 
model  constructed  from  full  waveform  tomography.  The  full  waveform  model  (Fichtner  et 
al.,  2009  and  2010)  is  derived  from  a  large  variety  of  observations  including  both  surface 
waves  and  body  waves.  They  also  implemented  a  realistic  3D  Q  model  by  Abdulah  (2007). 
Covellone  and  Savage  (2012)  also  used  a  3D  model  (Kustowski  et  al.,  2008)  computed 
from  a  combination  of  data  sets  that  included  surface  wave  phase  velocity  measurements, 
long  period  waveforms  and  body  wave  travel  times.  This  suggests  tomography  models 
derived  from  full  waveform  modeling  are  preferred  for  short  period  inversions. 

Since  our  results  indicate  the  use  of  3D  GFs  at  short  periods  has  limited  benefit  for 
the  particular  3D  model  that  was  employed,  a  more  attractive  option  to  evaluate  3D  models 
than  the  costly  3D  simulations  may  be  using  path-averaged  ID  velocity  models  derived 
from  the  3D  model.  We  compute  MTs  for  Little  Skull  Mountain,  HOYA  and  JUNCTION 
using  ID  averaged  velocity  models  for  each  source-receiver  path.  In  general,  the  averaged 
structures  have  a  Moho  ranging  from  30  to  37  km  and  a  very  thin  low  velocity  layer  (<1  km 
thick)  when  compared  to  the  ID  Song  et  al.  (1996)  model.  Here  we  present  results  for 
Little  Skull  Mountain  (Figure  45)  and  HOYA  (Figure  46).  Overall  the  MT  solutions  from 
20  to  50  seconds  Figure  45a-b  and  Figure  46a-b  are  similar  to  previous  ID  and  3D  MTs, 
the  differences  are  in  the  short  period  inversions.  For  JUNCTION  the  results  are  similar  to 
the  3D  inversion  except  the  path  to  BKS  cannot  be  modeled  by  the  averaged  ID  structure, 
and  the  solution  has  a  lower  ISO  of  52%.  For  Little  Skull  Mountain,  the  short  period 
inversion  with  averaged  ID  models  have  better  fits  to  the  data  at  ISA  and  CMB,  however 
the  resulting  mechanism  is  an  incorrect  oblique  strike-slip  earthquake  (Figure  45a).  The 
overall  lower  goodness  of  fit  with  averaged  ID  MTs  and  3D  MT  suggest  paths  from  Little 
Skull  Mountain  to  the  stations  need  further  refinement.  For  HOYA,  the  averaged  3D 
structure  is  actually  a  better  representation  for  the  paths  to  GSC,  ISA  and  PFO,  the  fits  to 
the  data  at  these  stations  are  higher  compare  to  the  3D  solution  due  to  better  agreement  in 
the  Love  waves;  whereas  again  the  paths  to  BKS  cannot  be  modeled  by  the  path-averaged 
ID  model.  ID  and  3D  comparisons  indicate  that,  in  many  cases,  well-calibrated  average 
ID  representations  of  the  Earth  structure  may  be  a  more  attractive  option  at  periods  as 
short  as  8  seconds,  but  3D  models  do  need  to  be  considered  as  shown  by  the  modeling 
results  for  the  path  to  BKS. 
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a.  2002  Little  Skull  Mountain  Moment  Tensor  Solutions 
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Figure  45:  2002  Little  Skull  Mountain  full  moment  tensor  inversion  results  with  path 
averaged  ID  Green’s  functions  at  different  frequency  bands,  (a)  Full  moment  tensor  focal 
mechanism  and  the  deviatoric  component  of  the  solution  are  plotted  as  well  as  the  moment 
magnitude  ( Mw ),  percent  double-couple  (DC),  percent  compensated  linear  vector  dipole 
(CLVD),  percent  isotropic  (ISO)  and  variance  reduction  (VR).  VR  from  deviatoric 
inversion  is  in  parentheses,  (b-c)  Data  (solid  line)  and  synthetic  waveforms  (dashed  line) 
plotted  from  left  to  the  right  are  the  tangential,  radial  and  vertical  components. 
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a.  HOYA  Moment  Tensor  Solutions 

FULL  MT  Deviatoric 

Composite  1 D  ( 

20- 50s 

Composite  1 D 
8- 20s 


b.  20-50S 


Mw 

DC 

CLVD 

ISO 

VR 

5.09 

23% 

26% 

51% 

83%  (82%) 

5.21 

10% 

42% 

48% 

45%  (42%) 

c.  8-20s 


Figure  46:  1TOY A  full  moment  tensor  inversion  results  with  path  averaged  ID  Green’s 
functions  at  different  frequency  bands,  (a)  Full  moment  tensor  focal  mechanism  and  the 
deviatoric  component  of  the  solution  are  plotted  as  well  as  the  moment  magnitude  ( Mw ), 
percent  double-couple  (DC),  percent  compensated  linear  vector  dipole  (CLVD),  percent 
isotropic  (ISO)  and  variance  reduction  (VR).  VR  from  deviatoric  inversion  is  in 
parentheses,  (b-c)  Data  (solid  line)  and  synthetic  waveforms  (dashed  line)  plotted  from  left 
to  the  right  are  the  tangential,  radial  and  vertical  components. 


3.4.7  Conclusion 

We  applied  source-receiver  reciprocity  to  compute  3D  GFs  using  the  FD  method. 
Using  the  full  waveform  MT  inversion  method  (Minson  and  Dreger,  2008)  we  analyze 
earthquakes  and  explosions  at  NTS  using  ID  and  3D  Earth  models.  Other  than  the 
computation  of  the  GFs,  we  applied  identical  data  processing  procedures  to  the  ID  and  3D 
waveform  inversion  and  use  the  same  station  configuration  to  allow  for  direct  comparisons 
between  the  source  properties  and  associated  uncertainties  of  the  two  models  and  evaluate 
the  results  in  different  frequency  bands.  Our  results  at  low  frequencies  show  good 
agreement  for  the  focal  mechanisms  between  the  two  models  and  slight  improvement  in 
waveform  fits  when  using  the  3D  model.  At  high  frequencies  the  advantage  of  the  3D 
model  is  limited,  mainly  due  to  poor  agreement  between  Love  wave  data  and  synthetics, 
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and  for  the  two  earthquakes  we  see  an  increase  in  non-DC  components  in  our  full  3D  MT 
results;  however  we  do  see  significant  improvements  in  the  3D  inversion  along  the  paths  to 
BKS  and  MHC  and  the  reduction  in  variance  is  especially  prominent  at  high  frequencies. 
In  addition,  we  see  better  agreement  between  Mw  and  the  reported  catalog  magnitudes 
when  3D  models  are  applied  in  the  inversion.  We  do  not  see  a  systematic  reduction  in 
uncertainties  associated  with  the  MT  when  3D  GFs  are  applied,  in  most  cases  the 
uncertainties  are  the  same  between  the  two  models  at  short  periods  but  vary  from  event  to 
event  at  long  periods.  The  3D  model  that  was  used  tends  to  add  a  false  isotropic  component 
to  MT  solutions  of  earthquakes.  While  this  is  a  negative  result  from  a  source-type 
discrimination  perspective  it  is  important  to  recognize  that  the  improvement  in  fit  afforded 
by  the  isotropic  component  is  not  statistically  significant.  Using  3D  velocity  models  for 
source  inversion  is  often  cited  as  a  means  to  improve  results,  however  this  analysis 
demonstrates  that  this  may  not  always  be  the  case  and  even  if  calibrated  3D  models  are 
employed  careful  analysis  of  uncertainty  and  solution  sensitivity  are  needed  before  non-DC 
components  of  MTs  of  earthquakes  can  be  interpreted.  The  results  also  show  for  the 
explosions  that  there  is  no  difference  in  the  ability  to  discriminate  the  explosions  from 
earthquakes,  which  is  a  positive  result,  although  reduced  uncertainty  in  the  source-type 
goodness  of  fit  space  using  the  3D  GFs  would  be  preferred.  The  results  however  do  show 
that  the  minor  non-isotropic  components  of  explosions  are  highly  variable  with  respect  to 
both  the  velocity  model  and  the  passband.  This  means  that  it  will  be  difficult  to  interpret 
such  results  for  mechanisms  of  non-isotropic  radiation  in  explosions  from  MT  inversions. 

Our  results  indicate  that  the  surface  wave  derived  3D  model  for  the  Western  U.S.  for 
seismic  MT  estimation  still  needs  further  refinement  along  the  paths  we  examined  (except 
perhaps  BKS/MHC)  to  model  wave  propagation  at  high  frequencies,  and  instead  of  the 
more  costly  3D  simulation  using  path-averaged  ID  models  for  short  period  inversion  is  a 
more  practical  option  in  many  cases.  A  3D  model  by  Shen  et  al.  (2013)  that  uses  additional 
constraints  from  receiver  functions  should  be  evaluated  and  compared  to  the  MT  inversion 
results  using  the  Moschetti  et  al.  (2010)  model  in  future  work.  In  this  study  we  have 
established  a  procedure  to  compare  and  evaluate  ID  and  3D  source  inversions.  When 
models  with  better  crustal  resolutions  become  available  we  can  begin  to  explore  depth 
sensitivity  for  both  explosions  and  shallow  earthquakes.  It  is  likely  that  when  trying  to 
improve  capabilities  for  other  regions  of  the  world  that  multiple  3D  velocity  models  will 
need  to  be  tested,  as  well  as  perturbations  to  those  models  in  order  to  evaluate  the  stability 
of  source  solutions  obtained  with  3D  velocity  models,  much  in  the  same  way  that  it  is 
necessary  to  do  so  with  ID  models. 

4.  CONCLUSIONS 

In  section  3.1,  we  define  elements  of  the  general  MT  in  terms  of:  (1)  its  nonnalized 
eigenvalues  that  characterize  its  source-type,  (2)  a  seismic  moment  scale  factor  that  scales  the 
nonnalized  eigenvalues  to  appropriate  size  or  scalar  moment,  and  (3)  its  eigenvectors  that 
specify  its  orientation.  We  utilize  this  formulation  to  implement  an  iterative  damped  LS 
inversion  scheme  to  invert  displacement  waveforms  for  best  fitting  eigenvectors  and  the  moment 
scale  factor  for  specific  eigenvalues,  which  results  in  the  best-fitting  MT  solution  for  a  specific 
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source-type.  However,  the  expressions  are  general  and  can  be  used  with  other  appropriate 
inversion  techniques  as  well.  Our  technique  is  successfully  demonstrated  by  estimating  best 
fitting  MT  solutions  for  synthetic  data  assuming  various  source  types:  cracks,  explosion  and  DC. 
For  low  frequency  displacement  waveforms  of  the  three  example  events,  we  find  the  NSS 
computed  using  the  inversion  approach  to  be  faster  and  more  accurate  by  VR  ~  0-3%  compared 
to  NSS  computed  by  forward  modeling  80  million  randomly  generated  MT.  To  better  constrain 
the  source-types  of  these  seismic  events,  we  employ  an  approximation  of  the  sign  function  in 
order  to  invert  FM  polarity  data  along  with  displacement  wavefonns  using  our  derivative-based 
inversion  scheme.  We  find  that  our  inversion  method  is  more  successful  than  the  random  search 
approach  in  recovering  the  MT  solution  with  the  maximum  VR  as  well  as  the  region  of  best- 
fitting  MT  solutions  or  source-types  for  NSS  computed  using  low-frequency  displacement 
waveforms  and  P-wave  FM  polarities  both  separately  and  jointly.  The  inclusion  of  P-wave  first- 
motion  observations  with  long-period  waveforms  narrows  the  range  of  possible  MT  solutions  in 
source-type  space  leading  to  improved  source-type  discrimination  (e.g.  Ford  et  ah,  2012;  Chiang 
et  al„  2014). 

In  section  3.2  the  effects  of  shallow  depth  of  burial  and  the  impact  of  vanishing  traction 
at  the  free-surface  on  the  recovery  of  the  seismic  moment  tensor  is  examined.  Theoretical 
vertical  dip-slip  Green’s  functions  associated  with  the  Mxz  and  Myz  components  have  vanishing 
amplitudes  at  shallow  depths  (<  1  km)  due  to  the  loss  of  traction  at  the  free  surface,  while  the 
Green’s  function  waveforms  look  similar  with  little  phase  distortion.  However  synthetic 
calculations  at  shallow  depths  show  moment  tensor  based  event  discrimination  is  reliable  and  the 
resolvability  of  the  moment  tensor  solution  depends  on  station  configuration,  noise  level,  the 
frequency  band,  and  the  velocity  model. 

We  are  able  to  recover  a  predominantly  explosive  source  mechanism  for  the  three 
HUMMING  ALBATROSS  chemical  explosions  from  moment  tensor  inversions.  Although  the 
source-type  uncertainty  analysis  shows  that  we  cannot  uniquely  characterize  the  events  as 
predominantly  explosive  using  only  wavefonn  data,  the  combined  wavefonn  and  first  motion 
method  enables  the  unique  discrimination  of  these  events.  This  method  has  been  applied  in 
previous  studies  which  have  shown  the  inclusion  of  P-wave  first  motions  in  addition  to  full 
waveform  data  eliminates  the  common  ISO-CLVD  tradeoff  (Ford  et  ah,  2012;  Chiang  et  ah, 
2014)  and  reduces  the  uncertainties  of  sparsely  recorded  underground  explosions  with  strong 
Love  waves  and  reversed  Rayleigh  waves  (Chiang  et  ah,  2014).  In  this  study  we  further 
demonstrate  that  incorporating  the  two  data  sets  is  particularly  useful  in  constraining  the 
isotropic  component  of  explosions,  and  the  method  not  only  applies  to  large  events,  but  also 
small  magnitude,  very  shallow  explosions  that  are  effectively  at  the  free  surface.  The 
combination  of  both  low  frequency  full  wavefonn  data  and  high  frequency  P-wave  polarities 
greatly  enhances  the  capabilities  of  the  moment  tensor  source-type  discrimination  method  in 
cases  of  sparse  station  coverage,  strong  Love  waves  and  free-surface  effects. 

The  moment  tensor  method  is  capable  of  event  discrimination,  and  although  yield 
estimation  using  the  recovered  absolute  seismic  moment  from  moment  tensor  inversion  remains 
challenging  and  can  have  large  uncertainties,  we  can  begin  to  put  error  bounds  on  our  moment 
estimates,  and  therefore  yield,  using  the  NSS  technique  and  combining  waveform  and  first 
motion. 
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In  section  3.3  numerical  experiments  were  conducted  to  investigate  the  effect  of  3D 
seismic  velocity  structure  on  the  recovery  of  the  seismic  moment  tensor  and  in  discriminating  the 
event  source-type.  The  results  of  the  synthetic  3D  model  sensitivity  tests  show  that  using  ID 
velocity  models  to  compute  Green’s  functions  for  use  in  long-period  (as  short  as  10  seconds) 
MT  inversions  is  possible,  and  that  at  least  synthetically,  we  do  not  find  significant  bias 
in  solutions,  nor  issues  in  being  able  to  discriminate  different  types  of  sources.  For  the  DC 
case,  we  can  recover  the  correct  mechanism  and  source  depth  using  a  universal  ID  velocity 
model  for  all  stations.  The  ISO  component  in  the  full  MT  is  less  than  20%  for  the  two  DC 
mechanisms  explored.  But  for  the  oblique  mechanism  the  inversion  recovers  a  solution  with 
a  highCLVD  component,  resulting  in  a  solution  that  is  50%  DC.  The  observation  of  false 
non-DC  components  up  to  50%  for  the  earthquake  source  shows  that  although  the  correct 
fault  plane  solutions  for  the  DC  can  be  recovered  at  relatively  long  periods,  the  CLVD 
component  can  be  large  in  full  MT  inversions,  however  the  NSS  remains  earthquake-like 
focused  at  the  origin  of  the  source-type  plot.  For  the  explosion  cases,  the  level  of  false  non- 
ISO  is  dominated  by  the  ISO-CLVD  tradeoff,  although  the  recovered  MT  consists  of  a 
CLVD  solution  the  DC  component  is  less  than  10%,  and  the  NSS  exhibits  the  typical 
explosion-like  signature  in  source-type  space  (e.g.  Ford  et  ah,  2010).  Including  additional 
data  from  first  motion  polarities  can  eliminate  the  trade-off  and  the  correct  solution  can  be 
recovered.  The  constrained  MT  solution  is  predominantly  explosive  and  has  a  very  low  DC 
of  9%.  Similarly,  results  for  the  explosive  composite  source  show  the  tradeoff  affects  the 
recoverability  of  the  MTs.  The  fault  plane  orientation  of  the  DC  component  in  a  composite 
source  cannot  be  constrained  well,  except  for  the  case  with  50%  DC  and  50%  ISO.  Results 
from  synthetic  tests  using  a  more  heterogeneous  3D  model  with  crustal  velocity  fluctuations 
following  the  von  Karman  distribution  show  the  full  MT  inversion  can  recover  the  correct 
DC  mechanism.  The  percentage  of  non-DC  remains  the  same  compared  to  the  result  for 
a  smooth  model  but  the  solution  has  higher  uncertainties  as  indicated  by  lower  waveform 
fits  and  larger  variability  in  the  NSS.  When  an  explosion  source  is  considered  the  false-DC 
component  increases  to  almost  50%  in  the  inversion.  The  NSS  exhibits  an  explosion-like 
signature  in  the  source-type  space  but  with  large  variability  that  extends  to  a  pure  DC  source. 
But  considering  the  best  solution  with  first  motion  constraints  shows  a  predominantly 
explosive  source  with  a  slightly  higher  DC  of  28%  compare  to  the  results  from  a  smooth 
model.  The  false  non-DC  components  in  earthquakes  and  the  false-DC  components  in 
explosions  arising  from  unaccounted  for  3D  path  effects  can  be  reduced  by  using  path- 
specific  ID  Green’s  functions.  The  Green’s  functions  are  computed  by  taking  the  average 
ID  model  for  each  source-receiver  path  from  the  3D  tomographic  model.  Path-calibration 
also  reduces  the  large  uncertainty  in  the  NSS  for  the  oblique  mechanism. 


In  section  3.4  source-receiver  reciprocity  is  used  to  compute  Green’s  functions 
accounting  for  3D  wave  propagation  effects  in  the  western  US.  These  3D  velocity  model 
Green’s  functions  are  then  used  to  examine  moment  tensor  solutions  for  NNSS  nuclear 
explosions  and  nearby  earthquakes.  Using  the  full  wavefonn  MT  inversion  method  (Minson 
and  Dreger,  2008)  we  analyze  earthquakes  and  explosions  at  NTS  using  ID  and  3D  Earth 
models.  Other  than  the  computation  of  the  GFs,  we  applied  identical  data  processing 
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procedures  to  the  ID  and  3D  wavefonn  inversion  and  use  the  same  station  configuration  to 
allow  for  direct  comparisons  between  the  source  properties  and  associated  uncertainties  of 
the  two  models  and  evaluate  the  results  in  different  frequency  bands.  Our  results  at  low 
frequencies  show  good  agreement  for  the  focal  mechanisms  between  the  two  models  and 
slight  improvement  in  waveform  fits  when  using  the  3D  model.  At  high  frequencies  the 
advantage  of  the  3D  model  is  limited,  mainly  due  to  poor  agreement  between  Love  wave 
data  and  synthetics,  and  for  the  two  earthquakes  we  see  an  increase  in  non-DC  components 
in  our  full  3D  MT  results;  however  we  do  see  significant  improvements  in  the  3D  inversion 
along  the  paths  to  BKS  and  MHC  and  the  reduction  in  variance  is  especially  prominent  at 
high  frequencies.  In  addition,  we  see  better  agreement  between  Mw  and  the  reported 
catalog  magnitudes  when  3D  models  are  applied  in  the  inversion.  We  do  not  see  a 
systematic  reduction  in  uncertainties  associated  with  the  MT  when  3D  GFs  are  applied,  in 
most  cases  the  uncertainties  are  the  same  between  the  two  models  at  short  periods  but  vary 
from  event  to  event  at  long  periods.  The  3D  model  that  was  used  tends  to  add  a  false 
isotropic  component  to  MT  solutions  of  earthquakes.  While  this  is  a  negative  result  from  a 
source-type  discrimination  perspective  it  is  important  to  recognize  that  the  improvement  in 
fit  afforded  by  the  isotropic  component  is  not  statistically  significant.  Using  3D  velocity 
models  for  source  inversion  is  often  cited  as  a  means  to  improve  results  however  this 
analysis  demonstrates  that  this  may  not  always  be  the  case,  and  even  if  calibrated  3D 
models  are  employed  careful  analysis  of  uncertainty  and  solution  sensitivity  are  needed 
before  non-DC  components  of  MTs  of  earthquakes  can  be  interpreted.  The  results  also 
show  for  the  explosions  that  there  is  no  difference  in  the  ability  to  discriminate  the 
explosions  from  earthquakes,  which  is  a  positive  result,  although  reduced  uncertainty  in 
the  source-type  goodness  of  fit  space  using  the  3D  GFs  would  be  preferred.  The  results 
however  do  show  that  the  minor  non-isotropic  components  of  explosions  is  highly  variable 
with  respect  to  both  the  velocity  model  and  the  passband.  This  means  that  it  will  be  dif 
ficult  to  interpret  such  results  for  the  mechanism  of  non-isotropic  radiation  in  explosions 
from  MT  inversions. 

Our  results  indicate  that  the  surface  wave  derived  3D  model  for  the  Western  U.S.  for 
seismic  MT  estimation  still  needs  further  refinement  along  the  paths  we  examined  (except 
perhaps  BKS/MHC)  to  model  wave  propagation  at  high  frequencies,  and  instead  of  the 
more  costly  3D  simulation  using  path-averaged  ID  models  for  short  period  inversion  is  a 
more  practical  option  in  many  cases.  A  3D  model  by  Shen  et  al.  (2013)  that  uses  additional 
constraints  from  receiver  functions  should  be  evaluated  and  compared  to  the  MT  inversion 
results  using  the  Moschetti  et  al.  (2010)  model  in  future  work.  In  this  study  we  have 
established  a  procedure  to  compare  and  evaluate  ID  and  3D  source  inversions,  when 
models  with  better  crustal  resolutions  become  available  we  can  begin  to  explore  depth 
sensitivity  for  both  explosions  and  shallow  earthquakes.  It  is  likely  that  when  trying  to 
improve  capabilities  for  other  regions  of  the  world  that  multiple  3D  velocity  models  will 
need  to  be  tested,  as  well  as  perturbations  to  those  models  in  order  to  evaluate  the  stability 
of  source  solutions  obtained  with  3D  velocity  models,  much  in  the  same  way  that  it  is 
necessary  to  do  with  ID  models. 
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APPENDIX 


Partial  derivatives  of  MT  elements  with  respect  to  (m0,  a1;  a2,  a3,  bt,  b2) 


For  partial  derivatives  of  a  quantity  y  with  respect  to  x  =  [a1}  a2,  a3,  blt  b2\,  we  will  follow  the 
convention 


dy  r  dy  dy  dy  dy  dy ' 
dx  .da1  '  da2  '  da3  ’  db1  '  db2. 


(A-l) 


Partial  derivatives  of  rx  and  r2 


drx  cti  drx 

— —  =  —  and  ——  —  0  for  i  —  1,2,3  and  j  =  1,2 
dat  rx  obj 


dr2 

—  [— ^ie23  >  ~^>2e23  >  (Pie21  +  b2e22)  ,  (a3e21  —  a1e23 )  ,  (a3e22  —  a2e23 )] 


(A-2) 


Partial  derivatives  of  ea 

ds-tj  8jj  a;  a; 

— —  = - —  for  i  —  1,2,3  and  j  =  1,2,3 

oaj  rx  rx 

deti  (A-3) 

— —  =  0  for  i  —  1,2,3  and  j  —  1,2 


Partial  derivatives  of  e2 


E  = 


F  = 


■  0 

0 

bi 

0 

0 

^2 

l-hl 

~b2 

0 

■  a3 

0 

0 

a3 

) 

—a1 

—a2. 

de2i 

daj 


e2i  ( dr2 \  Eij 

—  (  - —  )  H - for  i  =  1,2,3  and  j  =  1,2,3 

r2  \dajJ  r2 


de2i 

dbj 


e2i  ( dr2\  Fij 

—  (  — -  )  H - for  i  —  1,2,3  and  j  —  1,2 

r2  \  dbjJ  r2 


(A-4) 
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Partial  derivatives  of  e3  with  respect  to  [a1}  a2,  a3,  blt  b2]  can  be  computed  using  partial 
derivatives  of  ej  and  e2  by  simple  chain  rule.  For  example, 


e31  —  e12e23  —  e13e22 


de31 

dx 


de23  de12  de22 

—  ^12  ~ T  f  ^23  ~ T  ^13  ~ T  ^22 

dx  dx  dx 


de13 

dx 


(A-5) 


Partial  derivatives  of  with  respect  to  [a1,a2,a3,b1,b2]  can  be  computed  using  partial 


derivatives  of  e1(  e2  and  e3  with  respect  to  [av  a2>  a3,  bv  b2\  by  chain  rule. 


dMLJ 

dx 


dfi i 

dx 


+  eij 


den\ 
dx  ) 


+  d2 


(  9e2j  , 

y 21  ~ax  +  $2< 


+  ^3 


de3j 

dx 


+  e3  j 


for  i  —  1,2,3  and  j  =  1,2,3 


(A-6) 


The  partial  derivatives  with  respect  to  the  moment  scale  factor  are  straightforward. 


dMjj  2  M, : 

— - - - for  i  =  1,2,3  and  j  =  1,2,3 

dm0  m0 


(A-7) 
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